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1.0 SUMMARY 


Several aspects of a finite difference method for solving the unsteady transonic flow about 
harmonically oscillating wings are analyzed. The procedure is based on separating the 
velocity potential into steady and unsteady parts and linearizing the resulting unsteady 
differential equation for small disturbances. The differential equation for the unsteady 
velocity potential is linear with spatially varying coefficients. 

The analyses decribed herein concern methods of improving the accuracy and efficiency of 
the finite difference solution. The overstability of the current upwind differencing for 
supersonic flow is studied for the Klein-Gordon differential equation, which is the equation 
for the flat plate oscillating in supersonic flow. The operator is shown to be overly stable in 
that the finite difference solution is attenuated in the downstream direction exponentially 
in terms of the frequency and the grid size. A stable differencing is derived which has 
greater accuracy. 

The addition of a viscous term has little effect on extending the range of convergence of the 
relaxation procedure beyond the critical frequency. A simple downstream boundary con- 
dition is derived on the assumption that the vortex sheet dominates the flow on the down- 
stream boundary. The results obtained with this boundary condition are indistinguishable 
from those with the plane wave boundary condition. 

Difference equations are derived using an oblique coordinate system which aligns the 
coordinate lines with the leading and trailing edges of tapered swept wings. 

The additional terms required for convergence of row relaxation of three-dimensional mixed 
flow are also derived. 

An approximate method of aeroelastic analysis for high aspect ratio wings using a two- 
dimensional direct solution with a full three-dimensional steady-state potential is also 
described. 

Except for the addition of the viscous term and the revised downstream boundary conditions , 
the analyses presented here are yet to be implemented in the program for computing tran- 
sonic unsteady harmonic flow around airfoils. 



2.0 INTRODUCTION 


The purpose of the work described in this report is to continue the development of a means 
for calculating air forces for use in flutter analyses of three-dimensional lifting surfaces in 
the transonic flight regime. The work concentrates on a particular procedure which assumes 
small perturbations, the existence of a velocity potential, and simple harmonic motion, and 
uses finite difference theory to solve the resulting set of partial differential equations. The 
velocity potential is divided into steady and unsteady parts. The steady potential is calcu- 
lated using the classic nonlinear small perturbation differential equation. The unsteady 
potential is then calculated using a linear equation with spatially varying coefficients which 
depend on the steady flow. This study represents a direct extension of the research des- 
cribed in references 1 through 3. Reference 4 contains the latest results achieved in the in- 
vestigation while this report covers analyses which for the most part are yet to be imple- 
mented. The purpose of these analyses is to improve the efficiency and accuracy of the 
solution. 

Several different finite difference procedures are discussed in section 5.0. Subjects include a 
new operator for mesh points with supersonic flow which is stable but does not attenuate 
the initial data, the effects of adding a viscous term to the original differential equation on 
the convergence of the relaxation solution, and an alternative and relative simple down- 
stream boundary condition. 

Section 6.0 presents a method which uses a finite difference procedure over a limited inner 
region which is matched on the mesh boundary with an approximate linearized solution for 
the outer region. This has the two-fold purpose of reducing the number of points in the 
finite difference region and improving the exterior boundary conditions on the mesh. The 
derivation and a detailed set of equations are included in appendix D. 

Section 7.0 discusses two subjects related directly to three-dimensional flow. The first is 
an oblique coordinate system for swept and tapered wings. The second part discusses 
additional terms required to make row relaxation solutions converge when mixed flow is 
present. 

Section 8.0 discusses a finite span flutter analysis using the two-dimensional unsteady 
transonic program with a full three-dimensional steady velocity potential distribution. 
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3.0 ABBREVIATIONS AND SYMBOLS 


a Streamwise dimension of mesh region; also coordinate of downstream 

boundary 

b Root semichord of wing or semichord of airfoil; also vertical dimension of 

mesh region 

Cp Pressure coefficient, (p - p 0 ) / (1/2 P 0 U 0 "') where p is the local pressure, 

p Q the freestream static pressure, and p Q the freestream air density 

f(x,y,t) Instantaneous wing shape defined by z Q = 5f(x,y,t) 

f Q Undisturbed wing or airfoil shape 

f ] Unsteady contribution to wing or airfoil shape 

h Vertical mesh point spacing 

i.j,k x,y,z subscripts and indices for points in the mesh 

I,J,K 

i aFT 

k Horizontal mesh point spacing 

K Transonic parameter, (1 - M~) / (M“e) 

le, LE Leading edge 

M Freestream Mach number 


n,m Mesh point indices 

q / e - ico(7 - 1 )<p Q 

XX 

t Time in units of b / U Q ; also psuedo time defined by iterations in the complex 

differential equation for the unsteady potential 

TE, te Trailing edge 

U Q Freestream velocity 

x 0 ,y 0 > z 0 Physical coordinates, made dimensionless with the root semichord 
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x,y,z 


x ,y ,z 


x £e ,x te 


y 

P 


7 


AC 


P 


A<p l 


^‘te 

5 


e 



M 


$»*?>£ 


u (x,y) 




*1 

* 

CO 


Scaled coordinates (x 0 ,/i y 0 ,fxz Q ) for the three-dimensional problem; the 

scaled coordinates for the two-dimensional problem are x and y, with x being 
the direction of fluid flow 

Variables of integration 

Coordinates of leading and trailing edges 

T/Ky 

Vl -M 2 


Ratio of specific heats for air 

Jump in pressure coefficient 

Jump in at plane of wing or vortex wake 

Jump in ip j , at wing trailing edge 

Thickness ratio or measure of camber and angle of attack 
(5 / M) 2 / 3 
uM/(l -M 2 ) 

Critical value of Xj 

1/3 2/3 

Scale factor on y Q and z Q , n = 5 M 
Coordinates for swept and tapered wing 

Source distribution over mesh boundary for exterior panel method 

Complete, scaled perturbation velocity potential; also used for the unsteady 
potential in finite difference equations. With multiple subscripts, is used as un- 
steady grid point values of unsteady perturbation potential. 

Steady scaled perturbation velocity potential 
Unsteady scaled perturbation velocity potential 

Potential satisfying Klein-Gordon equation. With appropriate subscripts, 
represents boundary potential for matching inner and outer solution in appendix D. 

Angular reduced frequency (semichord times frequency in radians per second 
divided by the freestream velocity, cub / U) 
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4.0 FORMULATION AND SOLUTION 


A detailed mathematical derivation of the method for the solution of the unsteady velocity 
potential for the flow about a harmonically oscillating wing is presented in reference 1. The 
discussion here will be limited to a brief outline of the procedure for the two-dimensional 
flow. 

The complete nonlinear differential equation was simplified by assuming the flow to be a 
small perturbation from a uniform stream near the speed of sound. The resulting equation 
for unsteady flow is 

[k - (7 - 1 V t - (7 + iv x ]^ xx + ‘Pyy ~(2y> xt + </>tt) / e = 0 (1) 

where K = ( 1 - M ) / (mTc), M is the freestream Mach number of velocity U Q in the 

x -direction, x and y are made dimensionless to the semichord b of the airfoil and the time t 
to the ratio b/U Q . With the airfoil shape as a function of time defined by the relation 

y 0 = 8f(x,t) 

the linearized boundary condition becomes 

Y? y = f x (x,t) + f t (x,t) (2) 

The quantity 5 is associated with properties of the airfoil (such as maximum thickness ratio, 
camber, or maximum angle of attack) and is assumed to be small. The coordinate y is 
scaled to the dimensionless physical coordinate y Q according to 

y = S^M^yQ 

and e is given in terms of 5 by 

, 2/3 

e = (5 / M) 

The pressure coefficient is found from the relation 


C p = - 2e(^ x + <£ t ) 


The preceding differential equation is simplified by assuming harmonic motion and by 
assuming the velocity potential to be separable into a steady-state potential and a potential 
representing the unsteady effects. We write for a perturbation velocity potential 


'/> = ‘/>0 ( X y ) + ^l( x ’ y ) e 


icot 


( 3 ) 


and for the body shape 

y 0 = Sf(x,t) = 5 [f 0 (x) + f 1 (x)e 1C ° t J 


Since the steady-state terms must satisfy the boundary conditions and the differential 
equation in the absence of oscillations, we obtain 



(4) 


[K-(7+lV 0x ]*>o xx + v> 0yy '0 


with 


y>o y = f o< x > , y = o - 1 < x< i (5) 

On the assumption that the oscillations are small and products of may be neglected, 
equations (1) and (2) with the aid of equations (4) and (5) yield 

j[K-(7+ 1)^0 x ]^ 1 x | + ^l yy ‘ ( 2ioj / e )^l x + q ^l = 0 (6) 

where 

2 

q = oj / e - ioj (7 - 1 Vg 

XX 

subject to the wing boundary conditions 

¥>1 = fi x + kof](x) , y = 0 - 1 < x < 1 (7) 

A computer program for solving the steady-state transonic flow about lifting airfoils based 
on equations (4) and (5) was developed by Krupp and Murman (refs. 5 and 6). The output 
of this program or a similar program can be used in computing the coefficients for the 
differential equation of the unsteady potential. The similarity of the unsteady differential 
equation to the steady-state equation suggests that the method of column relaxation used by 
Murman and Krupp for the nonlinear steady-state problem should be an effective way to 
solve equation (6) for the unsteady potential Note that equation (6) is of mixed type, 

being elliptic or hyperbolic whenever equation (4) is elliptic or hyperbolic. Central differ- 
encing was used at all points for the y derivative and all subsonic or elliptic points for the x 
derivatives. Backward (or upstream) differences were used for the x derivatives at all hyper- 
bolic points. 


The boundary condition that the pressure be continuous across the wake from the trailing 
edge was found in terms of the jump in potential to be 


A n ~ A % e e 


-105 


( x-x te) 


( 8 ) 


where ^ is the jump in the potential at x = Xj. e just downstream of the trailing edge and 

is determined to satisfy the Kutta condition that the jump in pressure vanish at the trailing 
edge. The quantity At^i is also used in the difference formulation for the derivative <p i 

yy 

to satisfy continuity of normal flow across the trailing-edge wake. 


For the set of difference equations to be determinate, the boundary conditions on the outer 
edges of the mesh must be specified. In the original unsteady formulation, these boundary 
conditions were derived from asymptotic integral relations in a manner parallel to that used 
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by Klunker (ref. 7) for steady flow. A later formulation in reference 3 applies an outgoing 
plane wave boundary condition to the outer edges of the mesh. This boundary condition is 
numerically simpler to apply and, on the basis of limited experience, appears to provide 
equally good correlation. 

The preferred numerical approach to solving the resulting large order set of difference 
equations is a relaxation procedure, which permits the calculation to be made as a sequence 
of relatively small problems. However, as discussed in preceding NASA reports by the 
authors (refs. 2 and 3), a significant problem of convergence with the relaxation pro- 
cedure was encountered which severely limits the range of Mach number and reduced fre- 
quency for which solutions may be obtained. The authors currently feel the only practical 
technique for circumventing these instabilities is a full direct solution where the difference 
equations are solved “all at once” rather than by line relaxation. 
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5.0 ALTERNATE FINITE DIFFERENCE PROCEDURES 


Three different aspects of the finite difference formulation for two-dimensional unsteady 
transonic flow are examined in this section. The first concerns an analysis of the finite diff- 
erence operator currently applied to supersonic mesh points, the second examines the effect 
on relaxation solution convergence of adding a viscous term to the transonic equation, and 
the third investigates the application of an alternate downstream boundary condition. 

5.1 NINE-POINT OPERATOR FOR SUPERSONIC FLOW 


In solutions of difference equations it is important for the operators to lead to stable solu- 
tions. Zajac (ref. 8) has shown, however, that the usual upwind differencing for the wave 
equation leads to an overstable solution in which the solutions decay exponentially with x, 
the time-like variable. We have extended his results to apply to the Klein-Gordon equation 
The details are given in appendix A. For the flat plate in unsteady supersonic flow, the 
differential equation takes the form of a Klein-Gordon equation 


1 2 
V'xx- Y ^yy + x i ^ = 0 

where K = (m^ - l) / (m e), X| = o>M / (m - l), and e = (5 / M)^. 
The function \p is related to ipj by y?j = . 


(9) 


Equation (9) is the equivalent to the Helmholtz equation for subsonic flow and is seen to be 
hyperbolic. 


We assume that the region over which the solution is to be found is discretized by a uniform 
mesh in which the x spacing is k and the y spacing h, with the expression i p n denoting i// 

evaluated at x = nk and y = mh. Backward differencing on the second derivative with respect 
to x with central differencing of the y derivative yields for equation (9) 


2 2 2 
^n,m " 2 ^n-l,m + ^n-2,m " p (^n,m-l " 2 ^n,m + ^n,m+l) ~ k M ^n,m 

where p = k /(h l/K). In appendix A, an exact solution of the difference equation (10) is 
found in the form 


* 


imd^ ±inr n 


nm 


= e 


cos r 


) 


-If 2 2 2' 

where r = tan [_4p sin(d/2) i-k Xj 

given by 


" 1 1/2 

j . Similarly an exact solution of equation (9) is 


01 ) 


<// = exp(W m ) • exp(±ix n +\{ 


( 12 ) 


Note that this solution is oscillatory without damping. 

We compare equation (11) with equation (12) by setting Q = hy in equation (11), expanding 



( 13 ) 


in powers of h and k, and retaining only the first-order terms in h and k. This yields 


/ / 2 \ 

r / 2 \ i 

■ , . / v 2) 

k i ^ 2 

^nm = exp(mym) * expl±ix n ^— + Xj 1 • exp 

" 2\K +Xl / n 


We see that the difference equation solution has damping in the x direction and because of 
2 2 

the terms v and Xj the damping is greater for the higher frequency components in the 
solution and higher reduced frequency in the equation. 

A stable difference operator utilizing nine points instead of the usual five will eliminate this 
excessive damping. We shall use central differencing for the x derivative and for y we use 

— ^ a (^n+l,m+l " 2 ^ n+ 1 ,m + ^n+l,m-l ) + ^ ' 2a )(^n,m+l " 2 ^nm + ^n,m-l) 
h 

+ a (^n-1 ,m+l ' 2 ^n-l,m + ^n-l,m-l)J 


where a > 0 is the parameter to be determined to make the operator stable, 
leads to the following difference equation in place of (10): 

Equation (14) 

^nm " 2 ^ n- 1 ,m + ^n-2,m - p [ a (^n+l,m+l ~ 2 ^n+l,m + ^n+l,m-l 

) 05) 

+ (1 - 2a) ^ n rn +] - 2^ n)m + ’/ / n ) m-l) 
/ \ 

2 2 1 

+ a yl J n-\ ,m+l ~ 2 ^n-l,m + ^n-l,m-l/ “ 

k M ^nmj 

We assume a solution to equation (15) in the form 


i^nm - exp(infl) • expfima:) 

(16) 

In appendix B this is found to be stable for 


a> 1/4 

(17) 

and kXj < 2 

(18) 


Choosing a convenient value of the parameter a subject to equation (17) and k sufficiently 
small for a given reduced frequency and Mach number will thus ensure a stable operator. 
Expanding the equation resulting from equation (16) and retaining terms up to first order in 
h and k yields the solution 



We see that this is in exact agreement with equation (12). Note also that a influences the 
solution in the third-order terms in h and k and higher. A simple form of the difference 
equation results when a = 1/2, for then the middle terms in equation (14) vanish. The 
additional points will not affect the basic diagonal system either for the relaxation solution 
or the direct solution although computing the matrix is slightly more costly. 



5.2 ADDITION OF A VISCOUS TERM 

It was suggested that the addition of a viscous term to the differential equation might 
improve convergence of the relaxation process and extend the frequency limit for which 
solutions could be obtained. In subsonic regions the viscosity resulting from the first-order 
truncation terms in the difference equation does play a role in the convergence. Reinforcing 
this viscosity with an additional term seems like a logical approach to improving conver- 
gence. Accordingly, the following differential equation was investigated: 

^xxx + *^xx + ^yy + ^1 ^ ~ 0 
On the upstream boundary the conditions 

\p x = 0 = sin(7ry / b) 

were prescribed with = 0 on the other boundaries. The equation was differenced and the 
program coded along with the amplification factors obtained from a Von Neumann stability 
analysis. For the coefficient of viscosity ju set to zero, the relaxation converged for values of 
A I less than the critical values of 



where a and b are the horizontal and vertical dimensions of the mesh region, and diverged 

for Aj greater than Aj as predicted by the Von Neumann analysis. The additional viscosity 
c 

had very little effect on the convergence even when fairly large values of fx were tried. 

5.3 AN ALTERNATE DOWNSTREAM BOUNDARY CONDITION 


In the search for simple mesh boundary conditions which would improve the accuracy of the 
finite difference method for the unsteady subsonic flow over a two-dimensional flat plate, 
it was reasoned that, far downstream, the flow field is dominated by the vortex wake shed 
from the wing trailing edge. This boundary condition is easy to formulate since it depends 
upon the jump in potential at the trailing edge required to satisfy the Kutta condition. Thus 
each difference equation for the column of grid points next to the downstream boundary 
would contain four additional terms involving the four values of the potential in the neigh- 
borhood of the trailing edge upon which the jump in potential depends. The potential 
resulting from the wake is an infinite integral of a Hankel function and for two-dimensional 
flow is given by equation (109) in reference 1 in the form 


^1 ( X, Y 1 ) = 


where \|/ = exp^iA jM(x - x')J 



( 2 ) 

and H„ is the Hankel function of the second kind. 


+icoipj 


If instead of the pressure function 
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is prescribed on the downstream boundary, then the resulting integral in the equation 

0 

obtained from applying the operator — + ico can be integrated in closed form. From 

9x 

equation (C-4) in appendix C this is seen to be 

- iA ^te 


+ lCOipj 


exp^iXjM(x - x te )] -XjYK y Hj ( 2) ( Ajr ) / r 


2 2 

where r = Wx - x te ) + Ky and x^. g is the x coordinate of the trailing edge. 

The coefficients in the difference equations on the column adjacent to the downstream 
boundary for the potentials 


‘Pi,-! i -1 > ^1,-1 i ' ‘Pi, j -1 ' ‘Pi, j 

M 1,J m 1 M 1,J m M ,J m 1 


m 


where j is the y grid index in the row adjacent to the wing and wake and ij is the x grid 

index for the point at the trailing edge, are developed in appendix C. The equations were 
derived assuming that ipj is antisymmetric about the line y = 0, corresponding to the wing 

and wake, for the purpose of testing the concept as economically as possible. The resulting 
pressures on the wing differed insignificantly from the results obtained by assuming an out- 
going plane wave boundary condition on the downstream boundary. Furthermore, the 
pressures are not sensitive to the location of the downstream boundary. It therefore appears 
that the outgoing plane wave boundary conditions produce little reflection back to the 
airfoil although the distance from the wing to the downstream boundary is not so great that 
one should expect the disturbance to resemble very closely a true plane wave. 
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6.0 A FAR-FIELD MATCHING METHOD FOR TRANSONIC 
UNSTEADY FLOW USING THE DIRECT SOLUTION 


Chen, Dickson, and Rubbert (ref. 9) developed a method for matching the far-field 
boundary of the transonic steady finite difference mesh with an analytic outer solution. 
Their method has the advantage of imposing analytic boundary conditions at infinity while 
permitting a considerable reduction in the size of the mesh region. The mesh need extend 
outward only to where the flow is subsonic and linearized theory is valid, rather than to a 
distance at which the approximate evaluations of the outgoing wave boundary conditions, 
or alternatively the Klunker-type boundary conditions, are valid. The solution in reference 
9 was obtained using relaxation procedures. However, it is possible to obtain the far-field 
matching solution and the inner finite difference solution in a single step. 

Reducing the size of the solution will facilitate refining the mesh size and this is necessary 
for obtaining suitable accuracy with the finite difference method at higher frequencies. 
Also, the reduction in mesh points may be enough to make the direct solutions practical 
for three-dimensional problems. Alternatively, this matching procedure may provide better 
boundary conditions at the higher reduced frequencies, although the need for this is some- 
what reduced by the improved results presented in reference 4. 

The procedure is applicable to both two- and three-dimensional flow although the following 
derivation is for the two-dimensional problem only. In this section the basic ideas are 
sketched briefly but a detailed derivation is presented in appendix D. In section D.l, the 
basic integrals are discussed. In section D.2, the basic functions for the panel source distri- 
bution are presented and the form of the influence coefficient integrals defined. Far- and 
near-field approximations of the integrals are analyzed in sections D.3 and D.4. Subsequent 
sections formulate the boundary conditions, the wake, and the matrix coefficients in suffi- 
cient detail for coding into the direct two-dimensional solution. 

The derivation of the matrix elements are for a doubly symmetric grid distribution and 
symmetric steady flow so that the method may be evaluated as economically as possible. 
The basic integrals are simplified to the extent of requiring a single coded subroutine. In 
section D.l 1, the formulas for those matrix coefficients required by the outer solution are 
defined in a simple form suitable for coding. 

Following reference 9, an acoustic source distribution is prescribed on the outer edge of the 
mesh and a single vortex line imposed on the wake. The source strength is determined to 
satisfy continuity between inner and outer solutions of the normal component of the velo- 
city and the velocity potential at the outer mesh boundary. The vortex line accounts for 
the jump in potential of the wake. We assume that the velocity potential of the outer 
solution satisfies the linearized differential equation for the harmonic unsteady flow of a 
gas given by 

K<pi -2i(co/e)^i +ipi + (co / e)^i = 0 
‘xx X x yy 1 

From equation (109) of Ehlers (ref. 1), the solution to this equation given by a source 
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distribution on the boundary of the mesh and a doublet sheet from the wing trailing edge 
takes the form 

-a 114,. /* b 


iX< Mx 
& 1 


*1 


4i 


/ a iXj Mx r b 

[°u^u - °d^d] dx ' - ^ — J h - M'rjdy' 


-a 


( 2 ) 


+ Avji • x(x.y) 

1 te 

( 20 ) 


J ,2 _ 2 

(2) 

,/ 2 ,2 

Xj y(x-x) +K(y-b) 

h = H 0 

xi V( x - a i) + K (y - y ) 


where \ p u = Hq 

with similar expressions for \p^ and i p { . Here x is the potential induced by the trailing 


vortex sheet 


X = 

4i 


/ 


-iw(x'-l) , 

e \p v dx 


with \}J = H 0 (2 )( X] Y( X - x ') 2 + Ky" 


The quantity a is the source strength, a and b are the width and height of the mesh shown 
by the heavy line in figure 1, u and d denote upper and lower, and £ and r, left and right 
boundaries, respectively. For each 1 < i < i max designating the column for the upper and 

lower boundaries, and each 1 < j < j max on the side boundaries, we match boundary condi- 
tions on the normal derivative and on the potential at each boundary point with the finite 
difference solution. The number of values for the potential <pj from the finite difference 
equations is i max J max - 4. 


At each outer boundary point we assign a value of the source and construct a piecewise 
linear distribution of the source strength on the mesh boundary, utilizing for each boundary 
point the elementary singularity distribution in figure 2. The velocity potential in equation 
(20) after the integration takes the form 

‘max - * -imax"! 

V = £ Kn^un * °dn^dn] + £ [ ff 8n*£n ~ a rn^rn] + Al ^l tp ' X(x,y) ( 21 ) 
n = 2 n = 2 Le 

where the <p n terms are the functions of x and y resulting from integrating the basis function 

of figure 2 over the range x n _j to x n+1 or y n _j to y n+ j . The jump in the potential Ay> te 

is given by a linear combination of values of the potential y>| at points in the neighborhood 
of the trailing edge. 


We now match the solution of equation (20) with the inner finite difference solution. On 
the upper boundary we write for the velocity potential 




(u) 


(Vij +^ij _ 1 X/2 = F i (o,A<P\ ) 

V 1J max 1,J max V J \ He/ 


i = 2,3. 


Hnax"^ 


(u) 

where Fj (a,Acp te ) is a linear function of the o ' s and the \p's by equation (21) evaluated at 

the boundary point x = Xjy = b on the upper boundary. In the same manner, for the lower 
boundary, we obtain 



_Q|Q 



(d), 


(^il + m) / 2 " F i ( CT ’ A % e )’ i = 2,3,..., i max -l 
For the left boundary 

(£) 

(cp lj + ^ 2 j) /2 = F j ( a ’ A *l te )’ j = 2 ’ 3 ’-’ j max- 1 
For the right boundary 

(r) 

(^ i maxj + iPi max- 1 >j ) /2 = Fj C ^ te )’ J = 2 ’ 3 ’-’ W ' 1 

Similarly, we evaluate the normal derivative on each of the mesh boundaries and obtain: 


For the upper boundary 

Jmax '’ J max 


V\\ ~ *P[ j -1 

Umax ‘'Jmax 1 


*V = 


v\ - y, -i 

J max J max 


„(«)/ . V 

G i (“Ac, te ) 


' 2,3,..., i max -l 


For the lower boundary 

Vy = 


*i2“*il (d) . . ... . , 

Gj i 2,3,..., i max 1 


For the left side 

For the right side 
^x 


y>x 


y2-yi 


x 2 " X 1 


- G T( M *'*) 


, j 2,3,... ,j max -l 


j - i£>: i : 

'max- 1 ‘max" 1 ’- 1 _ (r) N . . . 

Gj J 2,3..., j max -l 


x i ' x i 
hnax max 


It is easily seen that the preceding systems of equations, along with the finite difference 
equations, yields 

'max jmax + 2 ('max + Jmax) “ ^ 2 

equations for the same number of variables to be determined. The increase in the number 
of variables required by the matching procedure is 


2 ('max + jmax) ~ ^ 2 

and is offset by the considerable reduction in the size of the mesh region. Because of the 
wavelike nature of the solution for the unsteady flow, using large mesh sizes near the outer 
boundaries of the finite difference mesh has been found to lead to poor representation of 
the flow field, resulting in inaccurate pressures on the wing. By decreasing the size of the 
mesh region, finer grids are possible with the same number of mesh points. The coefficients 
of the additional terms in the system of equations are derived in considerable detail in appen- 
dix D. 
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It seems worthwhile to make some assessment of the computer resources required to apply 
the inner and outer matching. In appendix D, the number of integrals was reduced by 
assuming a grid which has two lines of symmetry. For a grid with i max x points and j max Y 

points, the number of integrals Nj to be evaluated for the symmetrical problem is 

N I = [‘max + 2(im“ O ’ 2 ] 

Many of these integrals need not be calculated by actually performing the integration using 
the Bessel function routines. In regions where the grid is fine in the far field, the integrals 
for every second or third grid point need be computed with the intermediate points evaluated 
by interpolation. The integrals are more complicated than the coefficients in the potential 
solution and hence must be calculated more efficiently. Chen et al. (ref. 9) obtained a con- 
siderable reduction in computing cost as well as improvement in accuracy. Much of the 
reduction in cost will come from the smaller mesh region made possible by the better mesh 
boundary conditions. For the higher frequencies where the grid spacing must be fine even 
in the outer field, this smaller mesh region should result in a considerably smaller matrix 
equation to be solved. Unfortunately, the equations involving the values of the source at 
the boundary grid points contain nonzero coefficients for most of the source values. Hence, 
the use of the inner and outer matching procedure introduces to the matrix of coefficients 
a vertical strip of nonzero elements which nearly eliminates the banded property of the 
original finite difference matrix. 
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7.0 INVESTIGATIONS FOR THE THREE-DIMENSIONAL PROBLEM 


Three-dimensional investigations for this report were limited to two problems. The first 
concerned developing a coordinate transformation for swept wings that concentrated grid 
points in regions of large gradients of y>j . The results of this study, which are based on 
transformations used for steady state, are presented in section 7.1, and detailed derivations 
are presented in appendix E. Previous studies have resulted in a derivation of a coordinate 
transformation for swept but untapered wings (ref. 1) and development of a three-dimen- 
sional program using a cartesian coordinate grid (ref. 2). In reference 2 it was also shown 
that for the two-dimensional problem, row relaxation converged more rapidly than column 
relaxation but that additional terms were required for points at which the equation was 
hyperbolic in order for the relaxation to converge. These terms for the three-dimensional 
problem are discussed in section 7.2, with a detailed derivation presented in appendix F. 

7.1 AN OBLIQUE COORDINATE SYSTEM FOR SWEPT AND TAPERED WINGS 

The three-dimensional unsteady transonic flow program described in reference 2 utilizes a 
rectangular grid. Better accuracy with fewer grid points can be achieved by using an oblique 
coordinate system chosen to align the leading and trailing edges with coordinate lines and 
hence provide the capability of finer grid spacing along these edges. The transformation will 
also make the unsteady program more compatible with the steady program. 


In the same manner as Bailey and Ballhaus (ref. 10) we consider a transformation of the 
form 


x - xg e (y) 


a = 


c(y) 


= g(x,y) 


( 22 ) 


i? = y 


where c(y) is the chord of the wing at the station y and xg e (y) is the leading edge of the 

wing planform. Thus £ = 0 is the coordinate representing the wing leading edge while £ = 1 
is the trailing edge. 

The coordinates £,r? must be defined beyond the wing tip. To achieve this, the wing leading 
edge is extended all the way to the mesh boundary by a straight line having the same slope 
as the wing leading edge at the tip. To ensure that £,rj is single valued in the region beyond 
the wing tip, the trailing edge is continued analytically beyond the tip by a quadratic whose 
slope varies continuously from the trailing edge value at the tip to a straight line parallel to 
the leading edge extension, as shown in figure 3. Thus the functions xg e (y) and c(y) and 
their derivatives are defined over the entire r? range of the mesh region. 

Under the transformation of equation (22), the transonic unsteady differential equation was 
obtained in appendix E, and in conservation form can be expressed as 


17 




Outer mesh boundary 


Figure 3.— Oblique Coordinate System for Swept and Tapered Wings 



a 

dj 


(m“u + i + (v - 2ioo / eVj + ryj 


+ • 


9r? 


( 1 + C 7 c)y?i + iy> 

'r? (23) 


+ *Vl 


q + c' / c - (c' / c)' 


<Pl “ o 


where p = 1 / c(p) and v = - £c' / c - x'(rj) / c(i?) . A simpler nonconservation form is given by 

^ 3 3 

M — J7iu^, € - 2icoipj / e] + n — + + *l ff + q ^l = 0 (24) 

The conditions that the equation be hyperbolic for both forms is 


2 2 
H u + v <0 


(25) 


This is the condition that Bailey and Ballhaus used at first to determine when to employ 
upstream differencing in the derivatives. They found that upwind differencing for all super- 
sonic points was required to capture the shock. 


On the wing root plane we must apply the boundary condition of symmetry yq = 0. In 

l y 

terms of the oblique coordinates this condition becomes 

% + I ^l£ = 0 (26) 

This boundary condition is applied to the difference equation for points along the wing root 
and leads to some simplification. The boundary conditions on the wing and on the wake are 
unchanged under the transformation. 


Equations (23) and (24) may be differenced in the same way as described in reference 1 and 
formulas are presented in the appendix E. Because of the cross-derivative terms, the grid 
point pattern used to represent the difference operator contains the eleven points shown in 
figure 4 instead of the seven points for the operator in cartesian coordinates in figure 5. 


7.2 ROW RELAXATION FOR THREE-DIMENSIONAL FLOW 

In reference 2 it was found that row relaxation for the two-dimensional solution of the un- 
steady velocity potential converges more rapidly than column relaxation. When the flow is 
completely subsonic the same difference equations may be used for either row or column 
relaxation. However, for mixed flows, the row relaxation will diverge unless additional time- 
like terms are added at supersonic grid points. 


Following Jameson in references 1 1 and 12, we introduce the time-like variable associated 
with the iteration process in the form of 


(n) (n-1) 

^ijk "^ijk 


= At 



(27) 
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and obtain the following differential equation by taking the limit as At, Ax, Ay, Az go to 
zero in the difference equation about the point ijk (see eq. (F-8) in app. F). 


At 


Ay; 


Hx ' 2K " ' x + %y + ^zz ' 2 ^(%t + n zt + h*l t ] + W\ = 0 


By a transformation of t in the form 


(28) 

(29) 


r = ax + |3y + 7z + t 

equation (28) can be converted to 

K - 2i(w / eVj] + y? lyy + n zz ~ a jk*Vr + b jk^l T + W\ = 0 (30) 

where aj k > 0. Since u > 0 at supersonic points, the resulting differential equation is not 
strictly hyperbolic in r as it is for subsonic point. 


The terms yq and yq are truncation terms resulting from differencing the x, y, z deriva- 
r rr 

lives in the conventional manner. To render the yq term time-like we must add yq 

*xx *xt 

differences to change the sign of yq and a yq difference to cancel the yq term in equation 

l TT *t l T 

(30). The derivation is presented in detail in appendix F and the equation to be differenced 
is (eq. (F-19) of app. F) 


[uy> lx -2i(w/eV 1 ] x + ^ ]yy + ^ zz -2^ pi^i x +^l t ] + ^1 = 0 (31 > 


where 


h 


= 


1 + P jk > C > 1 


(32) 


h = (pi ( u x ■ 2i 60 / e ) / 2u ] - h 
$3 = ' Ay j [ b j + b k - ( a j + a k ) (r - 1 )] / r 


Ay j = y j+l " y j-l 

Pjk = Ay j / ( z k+l ' z k-l) 

r = the relaxation factor 

The factors in equation (32) for two-dimensional flow are obtained from equation (32) by 
dropping all terms with subscript k. 


Row relaxation has the same frequency limitation as column relaxation but its greater 
efficiency may make it worthwhile for frequency ranges in which it converges, while going 
to some form of direct solution for the higher frequencies. 
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8.0 FORMULATION OF AEROELASTIC ANALYSIS 


References 1 to 4 describe a practical procedure for calculating transonic air forces for 
harmonically oscillating airfoils. The frequency limitation problem discussed in references 2 
and 3 appears to have been overcome to the point where combinations of Mach number and 
reduced frequency of practical interest in flutter can be handled. The size capability of the 
pilot two-dimensional program has been increased to work with a practical number of mesh 
points for these analyses. Also, the solution program has been modified to treat multiple 
right-hand sides efficiently. However, due to the large size of the matrix inverse which is 
required, this procedure does not, as yet, appear to be practical for full three-dimensional 

configurations. The full three-dimensional problem involves an inverse of a 50 000^ to 

100 000^ T order complex matrix. This may eventually be practical through the use of the 
new vector machines, or through use of sparse matrix concepts. The following paragraphs 
discuss one use of the harmonic finite difference procedure in flutter analyses. 

It should be emphasized that the problem formulation provides superposable pressure distri- 
butions which can be used directly in conventional (e.g., V-g) flutter analyses. The flutter 
equations in matrix form and applicable to both two- and three-dimensional flows are pre- 
sented in detail in section 10.0 of reference 3. 

Use of the direct solution program of this report for practical two-dimensional flutter 
problems appears to be feasible. It is, of course, highly desirable to extend the harmonic 
analysis to full three-dimensional flow. However, a reasonable alternative may be to use the 
two-dimensional program to calculate the unsteady pressures at several spanwise stations 
with the equation coefficients being determined from the three-dimensional steady- 
state velocity potential. This would make use of the current capability and include the 
major three-dimensional effects of the shock and boundary layer through the steady-state 
potential, and could prove in the long run to be a valid economical alternative to the full 
three-dimensional calculation which would be much more expensive in terms of computer 
resources. The procedure may be summarized in the following steps: 

1 . Calculate the steady three-dimensional velocity potential distribution using a standard 
small perturbation program such as that of Ballhaus and Bailey. 

2. Use the two-dimensional unsteady program with the three-dimensional steady-state 
potential to calculate sectional harmonic pressure distributions at a set of spanwise 
stations. Using the steady potential ensures that the three-dimensional shock effects 
are incorporated in the results. 

3. Form a three-dimensional pressure distribution from the two-dimensional section dis- 
tributions. Additional finite span corrections could be introduced at this time. These 
corrections could be based, for example, on empirical data or steady-state analytical 
data. 

4. Calculate generalized air forces and perform flutter analysis. 



APPENDIX A 


OVERSTABILITY OF THE CANONICAL UPWIND SUPERSONIC 
OPERATOR APPLIED TO THE KLEIN-GORDON EQUATION 


A.l INTRODUCTION 

When partial differential equations are solved by numerical methods, an area of particular 
concern is the stability of the numerical operators employed. In the case of hyperbolic 
equations, in particular, it is required that the operators be stepwise stable, i.e., that errors 
at one stage are not magnified as the solution is stepped along in time (or in a time-like 
direction). Such stability may ordinarily be established by a Von Neumann analysis. 

It has been observed by Zajac in reference 8 that some operators may be so stable that the 
correct numerical solution is distorted by being attenuated in stepping along. He has called 
this phenomenon overstability. The situation here is that while the numerical solution will 
converge to the true solution as the step size is refined, for a given step size the error may 
compare poorly with that obtained using a less strongly stable operator. 

A. 2 ANALYSIS OF THE NUMERICAL SOLUTION OF THE KLEIN-GORDON 
EQUATION USING THE SUPERSONIC OPERATOR 

A.2.1 DEFINITION 


The Klein-Gordon equation 


* 


xx ~ k *^yy + Xl ^ ^ 


(A-l) 


where 


and 


K = (m 2 - l) /(m 2 c) , Xj =cuM/ 
e = (8 /M) 2 ^ 3 


(m 2 - l) 


bears the same relation to the flat plate equation for supersonic flow as the Helmholtz 
equation does for subsonic flow. Observe that when Aj = 0, the K-G equation becomes the 
wave equation as, analogously, the Helmholtz equation becomes Laplace’s equation. In the 
supersonic case, however, x and y are not treated identically in the discretization, but rather 
the time-like character of x is considered and a backward difference operator is used. 


A.2.2 DISCRETIZATION 


We suppose the region over which equation (A-l) is to be solved to be discretized by a mesh 
such that k is the spacing in the x direction and h is the spacing in the y direction. With the 
mesh point which is the n^ 1 in the x direction and m^ 1 in the y direction, there is associated 
a value i jj nm which is an approximation to i^(nk,mh), i.e., to the solution at this mesh point. 
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Using the backward and central difference operators in the x and y directions, respectively, 
equation (A-l) becomes 

2 1 2 2 

(^nm'“^n-l,m + ^n-2,m)^ j/(^n,m-l " “^nm + ^n,m+l ) I ^ ~^1 ^nni 


or. 


where 


7 7 

^nm- 2 Vl,in + V2,in = P"(^n,m-l " 2 ^nm + *n,m+l ) ~ k ~ X l *nm (A-2) 


p = k / (lfVK) 


A.2.3 STABILITY ANALYSIS 


The exact solution for the difference equation (A-2) is 

n im0 
*nm = “ e 

Substitution of (A-3) into (A-2) yields 

im0/ n n-1 n-2) 2 nf i(m-l)0 im0 i(m+])0"| 2 n im0 

e \a -2a + a /= p a Le - 2e +e J-k X.a e 


(A-3) 


From which on division by 


n-2 im0 
a e 


we have 


or 


(2 \ 2 2T 

(a - 2a + 1/ = p a I- 


2 9 2 

4 sin (6/2) - k~Xj 


] 


f 2 2 2 21 2 

Jj + 4p sin (0 /2) + k X j J a - 2a + 1 = 0 


(A-4) 


which is a quadratic equation in a. Since the coefficient of a is always > 1 , we can define 

7/2 


r 2 2 2 21 

cos r = [l + 4p sin (0/2) + k Xj J 

for 6 < t < it / 2 . Then (A-4) may be written as 

2 ( 2 \ 2 
a - 2\cos r/a + cos r = 0 

Solving for a we have that 


±ir 

a - e cos r 


(A-5) 


(A-6) 


(A-7) 


Thus since I a 1 < 1 for all r, the operator given in equation (A-2) is unconditionally step- 
wise stable. 

A.2.4 OVERSTABILITY ANALYSIS 

In this section we show that the difference scheme used in obtaining equation (A-2) from 
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equation (A-l) yields an overstable operator. 


First, from equations (A-3), (A-5), and (A-l) we observe that exact solutions of the difference 
equation (A-l) are given at mesh point (n,m) by 

n im0 


* 


nm 


a e 


or 


for any real 0 , where 


r = tan 


^nm e 
-1 


imd{ ±inr n ) 
\e cos t) 


(A-8) 


2 2 / 0 \ 22 

4p sin + k M 


y 2 


(A-9) 


Second, we observe that exact solutions of the differential equation (A-l) are given by 


*nm = ^(^m) ’ ex P\ ±ix n]/j + X 1 / 
for any real v. Observe that the solutions in equation (A- 10) oscillate without damping. 


(A- 10) 


Let us try to compare equations (A-8) and (A- 10). This is facilitated by letting 0 = hv in 
(A-8). For the first factor we have then 

im0 imhv 


= e 


= exp 


O^m) 


with y m = mh. This is the same as the first factor of (A-10). 


±inr 


Next consider the second factor of equation (A-8), e . Letting 0 = hv in equation (A-9), 
we have 


-1 


t = tan 


k 2(hv 
4 — ~ sin 
2 

h K 


, 2 2 
+ k X | 


X /2 


which for hv small yields 


-1 


r « tan 


2 2 

k v 2 2 

L “ir +k x i . 


! /2 


which for 


k small gives 


r « k 


v , 2 
/ i- Xi 
' K 1 


Thus e 


±inr 


exp^±ix n y— + Xj J , since nk = x n . This is the same as the second factor 


of equation (A-10). 

Finally, let us consider the factor cos r in equation (A-8). We have 


25 



cos r = exp(n log cos t) 


= exp<- -log 


( " 

exp<- - log 


2 2/hv 


2 ^ . 2 


1 + 4p sin I — l+K z Xj 




1 + k 
2 


2[ v 

3T 


. k [v ,2 
exp v nk • — l — + Xi 

2 V K 1 / 


• + Xi 


for small hi> 


for small k 


or. 


n 

cos r « exp 


k/ v_ 
2 \K 


+ X 


1 


'n 


Tims the solutions in equation (A-8) have damped oscillations. Note that the damping is 
greater for higher values of frequency v in the solution and higher values of reduced frequen- 
cy X| in the differential equation. 
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APPENDIX B 


A STABLE DIFFERENCING SCHEME FOR THE KLEIN-GORDON 
FORM OF THE FLAT PLATE EQUATION IN SUPERSONIC FLOW 

We now establish a difference scheme for equation (A- 1 ) and show that it is stable without 
introducing attenuation. As before, we suppose a uniform discretization in x and y such 
that x n = nk and y m = mh and denote the value of \jj at (x n , y m ) by \p nm 

The form of the difference equations to advance the solution from x n to x n+ j are obtained 
by the following substitutions: 


^xx( x n’ y m) ^(^n+l,m ^nm + ^n-1 ,m)^ 
^yy ( x n> y m)^ 2 j^ a n+ 1 ,m+ 1 - ‘^n+l,m + ^n+l,m-l) 


(B-l) 


(B-2) 


+ (1 -2a)^/ n rn+ j 2\p n m + +a (^n-l,m+l “^n-1 ,m + ^n-1 ,m-l)J 

where a is a parameter, a > 0, to be determined. Making these substitutions into equation 
(B-l) and multiplying by k^, we have the implicit difference equation: 

P ^(i/'n+l ,m+l _ ^^n+l,m + ^n+l,m-l)j| + ^ '^ a K^n,m+l ' ^nm + ^n,m-l) 

+ a (^n-l,m+l * ^^n-l,m + ^n-l,m-l) " C 1 ^nm " ^n+l,m " ^nm + ^n-l,m 

2k 2 222 

where p = — — and cj =k Xj 

Kh 2 

The parameter a is to be determined from a Von Neumann stability analysis, which we now 

ind imce 

perform. On substitution into equation (B-3) of y/ nm = e e and subsequent division 


(B-3) 


inf? imce 

by e e we obtain 

2f ice/ id -id) / -id 

p [ae \e - 2 + e /+ ( 1 - 2a) \e 

After using the identity 


J\ -ice/ id -id) ' 

2 + e /+ ae \e -2 + e / - 


-id) 


2 ice -ice 

C] = e -2 + e 


ix -ix 2 

e -2 + e = - 4 sin (x / 2) 


the preceding equation simplifies to 


1 r 


p [_- 4 sin (0/ 2)] [l - 4a sin (ce / 2)_ - c j 2 = - 4 sin 2 (ce / 2) 


Z 

Solving this equation for sin (a / 2) , we have 

2 2 2 
2 p sin (d / 2) + cj / 4 

sin (a / 2) = — 


(B-4) 


2 2 
1 + 4p a sin (d / 2) 

A necessary condition for stability is that equation (B-4) can be solved for real ce for every 
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real 0. This will he true if and only if 

2 2 2 
p sin (0 / 2) + cj / 4 

0< <1 (B-5) 

2 2 

1 + 4p a sin (0 / 2) 

The left-hand side inequality is automatically satisfied. The right-hand side inequality is 
equivalent to 


p~[sin 2 (0 / 2)] [1 - 4a] < 1 - cj 2 / 4 for all 0. 

If a is chosen such that a > (1 / 4), then the left-hand side is < 0 for all 0 (and p). Since 
the left-hand side is 0 for 0 = 0 we must have 

2 , 

1 - cj / 4 = 0 (B-6) 

2 

Using the definition of cj we have that equation (B-6) is equivalent to 


k\ { <2 (B-7) 

k 

Thus the difference equation (B-3) is stable for all p = — — - provided that 

(i) a > 'A hl/K* 

(ii) kX| <2 

Choosing a according to (i) for convenience, we can satisfy (ii) by selecting k sufficiently 
small for the given reduced frequency and Mach number. 


With these restrictions on k and a, we now find the solution ^ nm as h and k go to zero. 
Then, as before, we let 

0 = tu> 


and 


exp(im0) = exp(my m ) 


and equation (B-4) becomes 



Since h and k are small we have 
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a « 



akv 

1 + 

K 



The solution ^ nm to the difference equation then becomes 

^nm = e im0 e ina = exp(h-y m ) • exp^ix n ^-+X 1 ^+o(k 3 ) 

We see that this is the exact solution to the difference equation in equation (A- 10) and 
shows no attenuation of the initial value problem as the solution progresses through the 
mesh. Since a is the order of unity, its value affects only the third -order terms in the grid 
spacing. 


Before choosing a, it is convenient to write equation (B-3) in another form to maintain the 
generality. 


Tridiagonal Form 


Here we consider as known all terms whose i// superscripts are < n, and as unknown those 
terms with \p superscript equal to n + 1 . Thus equation (B-3) becomes 

“P a (^n+1 ,m+l“ ^n+l,m-l) + (* + ^ p ~ a )^n+l,m 
2 2 

= P (1 -2a)^ n+l m+1 + +P a (^ n -i,m+l + ^n-l,m-l) (B-8) 

+ [ 2(1 -p 2 (l -2a)) -cf]^ nm - (l +2ap 2 )^ n . ] m 

which represents a tridiagonal system for each fixed n. 
a = 1/2 


For this choice of a, equation (B-8) becomes 


2 (^n+1 ,m+l + ^n+1 ,n-l) + 0 + p )^n+l,m 
2 

y^n-Fm+l + ^n-l,m-l) + ( 2 ' C 1 )^nm'( 1 +p )^n-l,m 
For point relaxation this may be written 

^n+l,m ~ 2(1 + p) [^n-l,m-l + ^n-l,m+l + ^n+l,m-l + ^n+l,m+lj + 



^nm 


“ ^n-1 ,m 
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APPENDIX C 


EVALUATING THE WAKE INTEGRAL FOR THE 
DOWNSTREAM BOUNDARY CONDITION 


In the quest for a better formulation of the downstream boundary condition^ we assumed the 
unsteady perturbation potential on the downstream boundary plane to be dominated by the 
flow induced by the doublet sheet shed from the wing trailing edge. Hence, for two-dimen- 
sional flow, the velocity potential at a point (x,y) on the vertical downstream boundary is 

pivp.n hv oo 

exp[-iw(x'-x te )] • [</y(x-x', yi -yj)] _ Q dx' 

term of equation (109) in reference 1). The notation is that 

used in reference 1 . 




^te f 

x te 

(see, for example, the second 


Since \p has the form 
then 


i>~ ^(x-x',yj - y j) 

3i// d\p 

3y i 9yi 


and we obtain 


* 1 


/ oo 

exp[-ico(x'-x te )] • dx'j 


He 


where, from equation (1 13) of reference 1 
\ \)~ exp 

( 2 ) 

and Hq is the Hankel function. From (C-2) 

* 


[^iXjM(x-x')] • Hq^ Xj V(x - x') 2 + ^ 


yi 


] 


(C-1) 


(C-2) 


yi 


= j - exp [iX j M (x-x')] • X j y i H j fx j V(x - x') 2 + y^] J / V(x - x ') 2 + y j ' 


Since y j = Yk" y . we have 


( 2 ) 


yi 

where r ="\/(x - x') ^ + Ky^ 


i/'y =<-expiX 1 M(x-x') -XjYKyjHj (Xjr)Wr 


(C-3) 


In reference 1 on page 61, it was shown that the integral in equation (C-1) resulting from the 
combination ipj + it o</?] can be integrated in closed form. From equation (78) of reference 

( 1 ) this is seen to be 
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-i Ai/> 


tc 


+ 100iP| 


cxp^iXi M (x - x te )J • X j VK y j H 


( 2 ) 


(X,r) /r 


(C-4) 


where r = m - X(- e ) + Ky 


The jump in potential at the trailing edge can be found in terms of values of the perturba- 
tion potential at grid points in the neighborhood of the trailing edge. From equation (104) 
of reference 1, we have for points on the wing 

Ay, =*i Jm+ | -•#>i,j m -c s !^ i .j n ,+2-^i.j m+ l) 


- c. 


( (U) (L)\ 

^(^m’^ m -l)-( d sl F i + d s2 F i ) 


(U) (L) 

where the constants are given in equation (105) and Fj and Fj are the boundary con- 
ditions on the upper and lower side, respectively. Since for the sake of economy in comput- 
ing resources for the test we restricted our analysis to steady-state flows without lift, i is 
antisymmetric and 




m 


j +? - j _1 

‘Jni ,J m 1 


m 


(C-5) 


Then on the airfoil 


(C-6) 


A ^i = - Xj m + c sl(*ij m -l " ^i J m ) + c s2 (^i J m " ^i J m -l) + ( d sl F i (U) + d s2 F i (L) ) 

A ^i = -(c s l + c s2 + 2 Kj m + (c sl +c s2 )^j m .j + (d s iF i (U) + d s2 F i (L) ) 

where the constants are given in equation (105) of reference 1 . At the trailing edge the 
Kutta condition requires 

+ icoAi^j = 0 (C-7) 

at x = Xj j , from equation (37) of reference 1 , we have 

c 'i l ( Avi i,+r Av, ii l ) +d >i 1 ( A ' 1 ’'i 1 " A,, 'i r i) +iw ^i 1| =0 

where c i and d i are given on page 40 of reference 1 . Solving for Ac?; , t yields 
i i *1 1 

'(i " d ii, /c ii 1 ' i “/ c ii 1 ) + ( d i i , /c i il ) A ^i-i <c- 


8 ) 


Using equation (C-6) to define A^j ^ and A«^^_j yields 
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A *te : = [ ! ' (\ +i “) / C li J [- (c sl - c s2 + 2)^ iiJm 

, v / (U) (L)\ 

+ ( c sl + c s2)^i,j m -l " (^sl F i + ^s2 F ij J 


(C-8) 

+ ( d li, /c li,) '( c sl +C S2 + 2 K-lJm + ( Cs l +C ^)^l-l.Jm-l ’ ( d sl F i + d s2 F i f , L -1 
Hence we write 

^Pl tc = h 1'PiiJ m + >'2*1, j m -l +h 3VU m + h 4V,.i m -l + R <C ‘ 9) 


where 


(C-10) 


h l = -( c sl +c s2 + 2 )( 1 ~ d l ij l c l - -icj / C lij) 
h 2 = ( c sl +c s2)( 1 - d l ii /c l ij -ico/ c i. ^ 
h 3 = " ( c sl + c s2 + 2 )( d l il 7 C l il ) 


h 4 = ( c sl +c s2)( d lj i / c l i][ ) 

R =-(l -d,^ / c,^ -ia»/ci. ] ^d sl F ii <U) + d s 2 F ii <L) ) 

’ C 1 i ! ' ° 1 > 1 ) ( ds 1 F ‘ 1 - <iU> + ds2Fj 1 

We now apply equations (C-4), (C-9), and (C-10) as the boundary conditions on the down- 
stream boundary. Thus in difference form we write 

i -<A -1 j V 3 ) « +^j _i j 

‘max> J ‘max ‘max ,J ‘max 

+ 100 


x i ‘ x i -1 
max ‘max ‘ 


iAv 3 ! 


te 


exp 


[iX 1 M(a 2 -x te )]J-^— 1 H l (2) (Vj) 
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where a-? = 


( X *max + X i m ax-l) 1 2 md r J = V( X te ~ ^ + Ky J 2 • 


This has the same form as equations (1 19) of reference 1 with a simpler function replacing 

P- ; . Solving the previous equation for A<pi yields 
•max-* te 

¥>j. .(l+iw6 2 /2)-(l-iw5 2 /2V 1 . . = h A{ Pl F j (C-H 

hnax’l te 


where §9 = x ; - x= 1 
‘max ‘max 


Fj = -ijexp iX}M^a 2 - x te )J ‘XjT/KyjHj (Xj^j / (4rj) 


where 


%ax’J Ck3 %iax' 1 »i + C k4 A ^l te F j 

c k 3 “ (* ” ^^2 I 2 )/ (l + icu52 / 2 ) 
c k4 = 5 2 / ( 1 + iw5 2 I 2 ) 


Substituting for Ay?i yields 
*te 

V 7 } j ” c k3^i ~\ ] 
‘max’ J ‘max ‘ ,J 


+ c k4 Fj ^VPi v i m + h 2%j m -l + Vi r l J m + h 4Vl,j m -l + R 

In the difference equation for general i,j, the potential <P\. is replaced by the right- 

hnax - * 

hand side of equation (C-l 6 ) when i = i max -l , the x index of the downstream boundary 
plane. 
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APPENDIX D 


A PANEL METHOD FOR MATCHING THE OUTER SOLUTION 
WITH THE INNER FINITE DIFFERENCE SOLUTION 
FOR TWO-DIMENSIONAL TRANSONIC UNSTEADY FLOW 


D.l INTRODUCTION 


The velocity potential for the unsteady linearized harmonic flow produced by a source 
distribution on a line segment s is given by 

iXiMx . 

'\LI 


n 


4i 


-/• 


a(x')H 0 ' '(£)ds 


/ 2 2 

where f = Xi r , r = y(x - x') + K(y - y') . The derivatives take the form 


n 


2 

Xj K 
y 4i 


n 


x 4i 


jj/" o(x') [h/ 2 \d / f] (y - y') ds J e lXlMX 
J a(x') [ h j ^ 2 \f) / (x - x') ds je lXlMx + iXjM<pj 


For convenience, we shall introduce the cylinder functions 


(2). 


V“) =(r) V (?) 


2 ^i r ,2 i 

where u = (f / 2) = ~~ [^(x - x') + K(y - y' )J . The derivatives of the functions take 


(D-l) 


(D-2) 


(D-3) 


the simple form 


2 n(u) = _fi n+l(u) 


(D-4) 


and higher order derivatives are obtained by simple recursion formula derived from the 
differential equation. 

k+2 r 1 lc n 

C n (u) = -^(n + k + l)£ n (u) + £ n (u)J / u k>0 (D-5) 


To match the outer mesh boundary with the proper outgoing wave solution for the rectangu- 
lar mesh in figure 6, we prescribe the following source and doublet distribution. 
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*\ 


-iX |Mx _ 1 
~ 4i 


7 -£ yV()(u')tls 

/ [«u< x ' )V 0 ("») - ( \l (x X)( u d)] tlx ' 


I 

4i 


(D-6) 


~^j [»v<y'>v ( ,("j') o,( v')V (J )]‘ly' 

X | VK yA^ j ( j ^ 

8i y 


("o) dx ' 


X i,-H 


where a u , erg, ff r denote the source strength on the upper, lower, left, and right edges of 
the mesh region, respectively. The subscript i] + l denotes the point just downstream of the 
trailing edge. Accordingly, the u variables are defined by 


X] r 1 2 ] 

U u = “L( x -x')~ + K(y-b) J 


u', 


■[(X - x'f + 


K(y + b) 2 ] 


X, 


_a i) 2 + K (y -y') 2 ] 


(D-7) 


U i=- 


{(x- a 2 ) 2 + K (y- y ') 2 ] 


u b = ~T^ * *') 2 + Ky 2 ] 


To simplify the derivation we assume a symmetric configuration without lift in the steady 
flow. Then the perturbation potential for the unsteady flow is antisymmetric and 
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ag(y') = - CTg(-y') 


o r (y') = -o r (-y') 

Since we consider only the lower half plane we have 


(D-8) 


J f b /*0 

b °c(y') 2 o( u c) d y' = J_ h °£(y' ) [ £ o( u e ) - £ o ( U c + )J 


dy' 


(D-9) 


where 


± X 1 


U C = 


[( x -a i) 2 + K( y ± y') 2 J 


Similar relations hold for the right boundary with replacing a j . Finally, for <pj : 


-iXjMx i r a 2 




4i 


/ d 2 

°d(x')[ £ o( u d) - £ 0( u u)] dx ' 


(D-10) 


il{I h ff c(y')[ c o( u e + )- £ o( u c )] rf y' 


ro + 

J o r (y')[ £ o( u r + )- £ o( u r )] d y' 


^W e 


iXjMx 


where we have changed the sign on o r for convenience and ip w is the contribution from the 

doublet wake: oo . , r 

-iXjMx Xjf y V, 1+ i /• iw ( x ' x ij+l) 


V w ~~ 


8i 


7 


£ i( u o) dx ' 


(D-l 1) 


where 


x ij + l 


. A i r ,2 21 

u o = — k x - x > +K y J 

Note that equation (D-10) satisfies the requirement </>](y) = -</>j(-y) 

D.2 BASIS FUNCTION FOR THE FINITE ELEMENT METHOD 


For each station xj on the upper and lower boundaries we use a linear distribution of 
doublet strength with the source strength defined at these grid points. For the basis func- 
tion centered at x = x n , we follow Chen, Dickson, and Rubbert (ref. 9) and use 
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° I °n~ 


! ( x '' x n-l) 


n= (WM^ forx "-‘ <x ' <y " 


2 ( x n+l - x> ) 


for x„ < x' < 


where 5 k = x^ - x^j . This form was chosen so that 

f [a(x')/o n ]dx'=l 
n-1 

A similar relation holds for sources at station y n on the x = aj and x = a 2 boundaries. 

Instead of xi and x ■ we replace these values by 
'max 


x l+x 2 


and a 2 = 


x i + xj _! 
'max 'max 1 


respectively in the end basis functions. We treat the end points on the other boundaries 
the same way. 


iXjMx 


We now consider the perturbation velocity potential without the factor e . Then 

with 

iXjMx 

ipj = e \p (D-13 

we consider the contribution from the source distributions 

ip = \jj d + \jj% + \p r (D-l^ 

where \p a is the contribution from the integral over the lower boundary and and \}t are 
the contribution from the left and right boundaries, respectively. Substituting the basis 
function into the integrals and performing the integrations yield 

*max~* 

’/'d = E a dn^dn 

n=2 / n i c 


il'g- E »Cn^Cn 
n=2 


'Pr E a rn^rn 
n=2 


where 
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(D-16) 



where 5 ^ - y k - y ^.7 , k = 2,3,..., j m , and a similar relation holds for \jj m from the 
source distribution on the boundary x = ^ We note that, from equation (D-10), the 
potentials all have the form 


= ^(y) ~U- y) (D-iV) 

The integrals in equation (D-16) may be calculated by Simpson’s rule, requiring the evalua- 
tion of the £q(u) functions at five points for each integral. With efficient coding this re- 
quired evaluation of the cylinder functions at the mesh boundary points and at midpoints 
between them for each point the induced flow is to be calculated. Far-field and near-field 
expansions of the integrals also may be used to reduce computing costs. 

D.3 FORMULAE FOR A FAR-FIELD EXPANSION OF THE INTEGRALS 

When the distance from the center of the panel inducing the flow to the point x,y is large 
compared with the range of integration x n _j to x n+ j (or y n _j to y n +j) , then the functions 
£ n (u) may be approximated by an expansion of the form 

00 /I x 1 00 1 

£ n (u + Au) = X] ^ (u) • Au / k! = 2L, a^Au (D-18) 

k=0 k=0 

where u depends only upon x,y and points (x n ,bY (a, ,y n ) or (a 2 ,y n V We see immediate- 
ly that V 2 V /V ) 


a nO“V u) (D-19) 

a nl = c n ( u) = - Vl( u) 
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From equations (D-5) and (D-l 8) we see that the a nk satisfy the recursion relation 


[" a nk 

/[(k + 2)u] 


a n,k+2--[< n + k+1 ) a n,k+l + k+] 


or more conveniently 



T a n,k-2 

/(ku) 


a nk = -[(n + k-l) a n,k-l + 

(D-20) 

We need to evaluate the integration along x'. Thus we write 



fig(u') = Cg(u + Au) 


(D-21) 

°° k 

= L a 0k Au 




where 

a 00 = C 0 (u) 


a 01 =_c i( u ) 

a 0k = ” ” ^ a 0k-l + a 0k-2 / (k ~ 1)1 / (ku) (D-22) 


Substituting equation (D-21) into equation (D-l 6) yields 

1 


^dn o; 


E a 0k 


*(Vl +S „)m" UK |X-1 

/* x n+l k ) 

/ ( x n+l ~ x ') Au dx '/6 n +l > 
x n ) 


Au dx' / 5 n 


(D-23) 


where, as we shall see, Au = p (wj - x') (w2 - x') . For convenience we shall introduce 

the functions 


u n k (x) = Jx' n Au k dx' (D-24) 

Note that equation (D-23) contains only one of the functions in equation (D-l 6) for the 
sake of simplicity. 


We require the functions Ug k and Uj k , and for later considerations, For the lower 


boundary we choose 

u d' = u d + Au 


where 


u d 


2 

~~ [(x - x') 2 + K(y + b)j 


(D-25) 
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then 


Au= -~[(x-x') 2 -(x-x n )“] 
Au = ju [(w, - x'^W 2 - x')] 


where 


w i = x n , W 2 = 2x - x n , and /a = X j / 4 
Similarly, for the left-hand boundary we have 


ug = ug + Au 


where 


Un = 


u fi 


■[(x-aj) 2 + K(y-y') 2 ] 


(x-aj) +K(y 


-y,,) 2 ] 


from which 


Au = 


Xj K 


(w 1 -y')(w 2 -y')] 2 


Xj k 


then 


where wj = y n and w 2 = 2y - y n . If we define /a = ^ 

Au = m^(wj -y')(w 2 -y')] 

Au = /a [( Wl - x ’)( w 2 “ x )j 

and the integrals along the boundaries take the same form. 

We now evaluate the functions Uj^Cx') . Thus with equation (D-29) 

U 0k( x ’) -f ^ (w i - x ') k (w 2 - x ') k dx' 


Let wj - x' = | ; then 


U ok (x') = " M k / ^ ( w 2 ' w 1 + £) d £ 


Expanding the k td power of the term in parenthesis and integrating yields 


U 0k ( x ') = -M k E (j k )( w 2 " w l) k J ( W 1 - x ') 


, N k+j+l 


/ (k + j + 1 ) 


,'k\ J =0 

where I I are Newton’s binomial coefficients. 

Writing 


(D-26) 

(D-27) 

(D-28) 

(D-29) 

(D-30) 

(D-31) 

(D-32) 
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(D-33) 


/z^wj - x')(w 2 - Wj^ = ot 

( W 1 - x ')/(w 2 -wj) =t 


yields 


Now 


U ok (x')--a k ( Wl -x') E 

j=0 V ' 

U lk ( x ') = /*< /( w l -x') ( w 2 " x ) x ' dx ' 


t 


k+j + 1 


= w iu ok +M £( k ) 
j=o ’ 


k-j 


k k /k\ ( w 2 ~ w l) ( w l- x ') 


k+j+2 


k + j + 2 


(D-34) 


(D-35) 


Similarly 


k v 2 k /iA fJ 

= w,Uok+<* (w,-x')- E( j )inrjT3 

j=0 ' 


u 2 k ( x ') = -/ M k ^(w2- w l + ^) k ( w l -t) 2 d^ 


2 , k JS, /k\(«'2- w l) kJ ( w l- x ') k+,+2 

= w,U 0k + 2w,M E(J — 

k k /k\( w 2 " w i) ( W 1 " x ) 

rj\j/ k + j + 3 


(D-36) 


Note that wj - x = 0 for x = x n and Wj = w 2 for x = x n . For this special case the U nk take 
a simpler form 

r K / j \ ZK ; 

u okVM w r x ) dx 


k/ ,\ 2k+l , 

= - M (w j - x j / (2k + 1 ) 


U„ = 


k. 


. 2k+2 


lk = w l U 0k + M ( W 1 " x ) /(2k + 2) 


2 k , A 2k+3 , 

U 2k =_w l u 0k + 2w l^lk - ^ ( w l" x ) /(2k + 3) 


(D-37) 

(D-38) 

(D-39) 
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Denoting 


U mk(V = U mk 


(D-40) 


and substituting for into the right-hand side of equation (D-23) yields 
a 00 


1 w l 1 i 

4 . ~ v '* 


/ 5 r 


(D-41) 


f n+1 n+ll [ 

+ |_ x n+ 1 U 0k ' U lk J/ 6 n+l( a 0k 

Since the ag^ are functions of u, which we shall define for convenience in the form 

u = ji[(x-x n ) 2 + z 2 ] (D _ 42 ) 

and since the U nk are functions of ju, Wj = x n , W 2 = 2x - x n , we can define a general func- 
tion 

G 0 (M, x , z, x n _j, x n , x n+1 ) 


1 


a 00 _ y 

4i 2i ( S n+l + «n) ^ 


j f n-1 

n-1 “I 

|_ x n-l U 0k 

- U lk J 


/«, 


f n+1 n+ll ( 

+ L x n+l U 0k " U lk J /5 n+l( a 0k 


Then we have finally 


^dn = G 0VT ’ x ’ ^ (y + b) ’ x n-l’ x n’ x n+l; 


, 2 

P'1 K 


^£n Gq\ a > y> ( x ' a l) / yn-1’ ^n’ ^n+l; 


^rn “ G 0' 


4 

/ 2 
M K 


> y> ( x -a 2 )/tK,y n -Dyn> V n +L 


(D-43) 


(D-44) 


D.4 FORMULAE FOR THE NEAR-FIELD EXPANSION OF THE INTEGRALS 


For the near field, the argument of u in the £ n (u) functions is assumed sufficiently small 
that the power series of the functions may be integrated term by term. Now 


J n «) = 



, 2k 

a/ 2) (-D 


k 


k!(n + k)! 
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(D-45) 


k=0 

n-1 2 

where i//(l ) = 7 and t|/ n = - 7 + 8/m Since u = (f / 2) then 

m=l 

„ k k „ 

8 o w = J otf ) = E — r = £ c ok u 


Similarly 


k=0 (k!) k=0 

or, k k „ 

«1 (u) = 2j,(0/f= E 77777 = E c ik u 

k=0 k=0 


k 2 

where c ok = (-l) / (k!) and c lk = c Qk / (k + 1) 


We also have 


Y 1 J 2 00 k 

8 0 (u) = Y Q (f) = - log u 8 q (u) - - £ lK k + 1 ) c 0k u 

7 ‘ 77 k=0 

8 0 Y (u) = Y 0 a) = - log u 8 0 J (u) - - V k + l)c ok u 

ir tt 

k=0 


8, Y (u) = 2Y itf) /? = -— + — log u 8j J (u) 

1 1 TT1 1 TT 1 


7TU 7T 


- ~ E l<K k + l) + ^(K + 2)fc lk u 

77 k=0 


Since 8q(u) = Hq^V) = 8q J (u) - i8Q Y (u), then 


00 \ lc 0k I k 

W u)=: E k Qk - logu[u 

k=0 ' 


Tt 


1 + — i//(k+ 1) 
7T 


(D-46) 

(D-47) 

(D-48) 

(D-49) 

(D-50) 

(D-51) 


(D-52) 


where 


b 0k " c 0k 


(D-53) 



Similarly 


where 


«i(u) = E 


k=0 


1C 


lk 


>1k lo 8 u 

7T 


k i 
u + — 

7TU 


b lk = c lk 1 + ~[^ k+ D+ 0( k + 2)], 
We need to evaluate integrals of the form 

/ ,n k 
x u dx 


(D-54) 


(D-55) 


and 


w nk ( x ') - J' x' u log u dx' 


(D-56) 


For integration along x', we have 2 


•[(x-x') 2 + K(y + b) 2 ] 


We now factor the quadratic in a form similar to equation (D-29). This leads to 

2 

X, 


u = 


(w 0 -x')(w 0 -x') 


where wq = x + i(y - b)T/K and wq is the complex conjugate. Similarly for integration 
with respect to y' we have 


2 

Xj K 


where 


As before we define 


[(wj -y')(w r y')] 


w 


j =y + i(x-aj) /1/K . 


u = ju[(w - x )(w - x )] 


(D-57) 


and obtain a general formula that holds, with appropriate arguments, for each of the three 
boundaries. Thus 


V, 


Ok 


. r k H , 

- I u dx 

= pk J (w - x') (w - x') k dx' 


Comparing with equations (D-34) and (D-36) we obtain 

v k ^ /k\(w - w) k J (w - x') k ] 1 

v ok ~ ~ v 22 (j 

j=o 


k+j + 1 


Similarly from equation (D-35) and (D-36), we have 


k 


V lk = wV 0k + ^ E 


j=0 XJ 


k\(w-w) (w-x') 


k-j , k+j+2 


k + j + 2 


(D-58) 


(D-59) 


(D-60) 
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V 2k ~ 2wV lk 



(- ^k-j ,,k+j+3 
(w - w) (w - X ) 

k + j + 3 


(D-61) 


Now 

then we have 


log u = log n + log(w - x') + log(w - x') 


u log u dx' 

/*V ok (x')dy' _ _ /*Vok(x')dy' 

= V 0k logM + V 0k log(w-x)+y — w _ x — + V Qk log(w - x ) +y - - 3 "_ x> ■ 

Substituting for VQ k yields 


W 0k = V 0k lo 8 M + 2 Real< V ok log(w - x') 


k / . _ k-j , k+j+1 
k —, / k\(w - w) (w - x ) 

+ M j! ~— 


j=0 VJ ' 


(k+j + 1) 


Now 


w 


/ J^ 

x'u log u dx' 

/•V lk (x')dy' 

= V lk logju + V lk log(w - x') + J • 


w - x 


Vi k (x')dx' 

+ V lk lo ^ - x ') V w-x' 


Substituting for Vg k yields 


w lk( x ') = Vj k log n + 2 Real j V lk log(w - x') 

u / \ - k-j , k+j+1 
k _ / k\(w - w) (w - x ) 


(k + j + 1) 




i k /iw— >3, ,.k+j+2 

k / k\(w - w) (w - x ) 


-<* E L 


j=o 


J 


(k + j + 2) 


(D-62) 


(D-63) 
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Similarly 


ial < ’ 


w 2k( x ') = v 2k ] °8 M + 2 Real j V 2k log(w - x') 

2 k /k\ (w - w) k \w - x') k ^ 1 


+ w M E ( , 

j=0 


(k + j + 1) 


(D-64) 


j=0' j ' (k+j + 2) 2 


i k / \ - k-j k+j+3 

^k^/k\(w-w) (w-x') 

j=0' J ' (k+j + 3) 2 


With the functions Vj k and Wj k defined and with 

V jk( x n) = Vn jk 

Substituting equation (D-64) into (D-23) leads to 


*dn “ ■ 


+ 


^dr 



-( v ik 


1 n\ 

- V lk /S n + 


/ n+1 n\ / n+1 n\ 

x n+l \ V 0k " V 0k / - \ V lk v lkj 


I 5 n+l 


OO 

E c 


2 *( 6 n+l + 8 n) k=0 


Ok 


/ n-1 n\ 

x n-l \ w 0k " w 0k / 


/ n-1 

n\l 

- V w lk 

- w lk^J 

/ n+1 

n\l 

- V w lk 

” w lk/ 


f n+1 n \ 
x n+l \ w 0k " w 0k ) 



Writing for the variable w 
and defining the function 


w = x + iz 


P 0 (m> x > z > x n-l’ x n> x n+l) 


(D-66) 


as the left-hand side of equation (D-65), we obtain 

L 2 

^dn “ P 0y 4 ,x,1^(y + b),x n _ 1 ,x n ,x n+1 y (D-67) 

L\ 

^n = p 0\ v “^~' »y>( x -ai) / Vk, y n -i,y n »y n +i^ 

and similar relations for i// un and \p m . 

D.5 CALCULATION OF NORMAL DERIVATIVE TO MESH BOUNDARIES 

The calculation of the contribution to <pi on y = -b from the source distribution on the 

y 

line y = -b is very simple. We consider 

■ a 


1 r 2 

t//(x,y)=-— / a(x')%(u)dx' 

2 41 at 2 

Af r 2 2l «— A i f 2 2~1 

where u = |(x - x') + K(y + b) J . Let IfK (y + b) = 77 then u = ~^~[(x - x') + 77 J 

and = "Vk xJj 1t j 


*'l V /»“2 

J a(x')«j(u)dx' 


d l 

Near u = 0, fij(u) = i / 7ru + 0(1) ; thus near 77 = 0 the integral takes the form 

• a 2 


1 - 77 

^ 2tt 


/ 


a(x )dx 


, 2 2 
1 (X - X ) +7} 


Since the integral does not exist at 77 = 0, we introduce the variable 

(x' - x) / t? = | 

,(a 2 - x ) / 77 


and obtain 


^ 2tt 


/ 




a(x + Tj£)d£ / 

“'(aj-x) / 7? 

If y goes to -b through values of y > -b, then 7? is positive and for aj < x < a 2 , the limit as 
goes to 0 becomes 
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and 


'P r) = o(x) I 2 


\jj y = VKa(x) / 2 


If y goes to -b through values of y < -b, then the limit becomes 

^y =- 

Similarly for sources on the y' axis we write 


* y = -VKa(x)/2 


- 


^2 

TT f cr(y')C 0 (u)dy' 
i b i 


where u 


'1 


(x-aj) +K(y-y') , let ifK (y' - y) = i?', then 




■L 


■lK(b 2 -y) 


4iVK^ CT(T? ') e o(u) d ^ 


(D-68) 


(D-69) 


1/K(b r y) 

Since this has the same form as the a(x') contribution, we see that 

* X = <K y)/(2VK) (P-70) 

for -b < y < b and x going to aj through values of x > aj . For x going to aj through values 
of x < aj, we obtain 


IK = - a(y) /CZVK) 


(D-71) 


D.6 THE BOUNDARY CONDITIONS ON THE MESH BOUNDARIES 


To match the interior finite difference solution with the outer finite element solution, we 
make the values of the potential from the two solutions and the values of the normal deriva- 
tives from the two solutions equal on the mesh boundary. Thus on y = -b = (y j + Y 2 ) / 2, 
we match the values of and 


from the two solutions while on x = aj 
we match the values of 


(x 1 +X9)/2 and on x = a-, = ( x= + x- i \ / 2 , 
\ 1 z / 1 \ hnax J max / 


from the two solutions. 


and i^j x 


We could actually evaluate v? x and «/> y for the outer solution by differentiating i p from the 
wake and source distributions; but to simplify the programming, we will approximate 
at x = Xj by + ^ 2 ) / 2 on the lower boundary both for the finite difference solution 
and for the finite element solution. For the derivative with respect to y, we take 

^i y = fc'^ii) 1 (y 2 ~y\) 
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for the finite difference solution and from the exterior solutions for the source distributions 
on the sides x = ai and x = ao ; but we will use 

limit ^d 

y -*■ - b ^y 

from the source distribution terms for the lower boundary y = -b. 


Let y?£, i p r , i be the contributions to the exterior solution from the left, right, lower 

boundary source distributions, and the wake. Then the boundary conditions on y = -b at 


x = Xj are 


+ + 

*i 1 + ^i2 1 + ^fii2 + ^ri 1 + ^ri2 + ^di 1 + *>di2 


' ^dil + ^di2 *wil + ^wi2 

- + - 


(D-72) 


2 2 

^il'^i2 *fiil -^£i2 + ^ril “^ri2 


yi-Y2 




+ lim 
y - b 


9y 


^dil " ^di2 Avil ” ^wi2 
+ 


V1-V2 Vl -V2 

where the + and - denote the lower source distribution and its image at y = b, and 
+ + 

y? r - <p r - v? r and - ¥>$> . Since we are interested in the outer solution at y = ~b 

through values of y < -b, then from the preceding equation we obtain 

*il -^2 = ^Pgil ~ *gj2 + Irfl ~ ^ri2 
2 


(D-73) 


iXjMxj 


-VKa di e (y i - y 2 ) / 4 - 

Adding and subtracting equations (D-72) and (D-73) yields 


*dil “^di2 


/ + + v iX i Mx i (yi-y 2 )^ 4+ ^wn 

m - nw + *ril - *dil + V^dil + y>di2 ) / 2 - a di e (D-74) 


■ / + + V iXjMxi (yi-y 2 )W/4 + ^ wi 2 

^i2 - + ^ri2 ■ ^di2 + V^dil + ^di2/ / 2 + ff di e 


(D-75) 


iX^Mx 


Because of the factor e , the normal boundary conditions on x = aj and x = ^ take a 
different form. For the derivative with respect to x we have 

+ 

iXjMxr + + -i 

— =e + iXjM>//£ J 
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Since we must approach x = aj from values x < aj, we obtain from equation (D-71) 

iX i Ma i 

/ 2YK 


lim + 

iXjM^g -ag(y)e 


x^-aj 3x 

Thus the boundary conditions on x = aj become 

+ + - 

^ij + v?2j n 1 j + ^£2j - n 1 j - ^£2j + ^r2j + ^r 1 j 


+ 


^dlj + ^d2j ,^wlj + ^w2j 


(D-76) 


^lj~^2j = - ^glj + mj + ^rlj~^r2j ^wlj'^w2j 


xi-x 2 


X 1 " x 2 
^dlj " ^d2j 


X 1 ' x 2 


iX jMaj 


(D-77) 


- °8 (yj) e ’ " 1 /dV k)+ iX[M ( n ' J + ~ ) 


x, -x 2 

Adding and subtracting equations (D-76) and ( x j - x 2 ) / 2 times equation (D-77) yield 
/ + + \ - / iXiMai \ 

^lj-<W£lj +^2j) -^lj + ^rlj + ^dlj + ^wlj-^ je 1 /(4 Yk|(x 1 -x 2 ) (D . 78) 

/ + + \ - ( iXjMaj 

*2j = “o(*81j + ^2j) '^>2j + ^r2j + ^d2j + ^w2j + \° r £j e /(4 Vk^/ (^x j - x 2 ) (D-79) 

where a 0 = 1 / 2 + iMXj(xj - x 2 ) / 4. 

On the right-hand boundary we apply the boundary conditions 

yf\ j ~ Vi _] j ) V-’JJi i - V^£i -1 j d j - V^(Ji -J i 

\ 'max- 1 'max 1>J / X 'max j X 'max u 'max J u 'max '> J 


( 


x i - x j 
'max 'max 


0 


( x i 


- X; 


i 'M -1 

'max 'max 


) 


(D-80) 


V^ri i " ^ r i _j j ^vvi i _ *^wi -1 j i ■ 9^5 

“max J "max 1,J w 'max J w 'max ' ,J lim u ^r 


+ 


max 'max 


max J 'max 


max J 


max 


+ + 

°r\ i ^ri 
"max J "max 


-1 x i ~ x i -1 

'max 'max 

x^ a 2 

3x 

I i ^di i ^di - 1 i 

'max J u 'max 1,J 

2 


(D-8 1 ) 

- 1 i ^ri i ^ ^ri - 1 i 

; 1,J max J "max 1,J 

^wi i 

WA max J 

^ w 'max" 


Now 


+ 

lim ^r 


x ^ a 2 3x lX l M ^r + x 


+ lim iXjMa 2 3 <p v 


a 2 


3r 
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(D-82) 


Since we approach a 2 from the outside, or x > a 2 , we have from equation (D-70) 

lim 9 *r + iX i Ma 2 y ^ 


lim 0 ^r + iX l Ma 2 ( _v 

-»■ a 2 “ = i^i^r + e a rj ^Yk) 


x ^ a 2 9x " 1 « A ' 

Thus equation (D-80) with the aid of equation (D-82) becomes 

P\ j - 'Pi -I j Pjli j ” P%\ -1 i *^(Ji i ~ Pd\ -1 i 

‘max- 1 ‘max ‘’- 1 Xi max J *‘max ‘> J UJ max J “‘max ‘’- 1 

x i ~ x i -1 x i ~ x i -1 

max max 1 ‘max ‘max ‘ 


Pri i - Pri 

“max J “i 


l max 
‘max'l’ J 


(D-83) 


x i x i - 1 

‘max ‘max 


T + 

Pri i Pri -1 i 
„ max J “max ‘ ,J 
+ iXiM 

i 9 


iXiMao ^wi i -l ^wi -1 i 

+ e a rJ - /(2 vK)+ 

and the other boundary condition is 

P\ i + P\ -1 i ^fii i -1 i + ^di i *^di -1 i 

‘max J ‘max ‘ > 3 ‘max J “‘max u ‘max J u ‘max ‘ ,J 

2 2 

^ J_ t- 

'rL ; + V” 


(D-84) 


+ + + 

™ ^ri i “ *Pr\ -1 j ^ ^ri i - 1 , 1 * 

Ai max J AA max AjJ AA max J AA max ,J 

H — 


^ w ‘maxj ^ w 'max ^ 5 j 


Adding /x; - x ; _i\ / 2 times equation (D-83) to equation (D-84) yields 

y ‘max ‘max f ' 

*P\ -1 i “ -1 i + ^di -1 i " V? r i -1 i + ^wi ”1 i 

A max AJ * A max A?J UA max A > J AA max AjJ WA max A ’ J 


+ + \ iXjMa 2 

+ 0i l(^ri i + V^ri -1 i~ e a - ,x 

M A1 maxJ A1 rrmx 1 >J/ 


(D-85) 


rii-i x i -1 1 /(4 Yk ) 
max’ “max 1,J / iJ V ‘max ‘max 1 ' 


Subtracting leads to 


Pi \-p£i i + PA\ i - Pri i + ^wi i + a l(^ri i + Pr\ -1 i 
‘max J max’ u max J max J W1 ma x J 1 i “maxJ “max 


max J 


'max’ “max ’ 


+ e 


iXjMa 2 


<i A /(4VK) 

ox/ 1 I 


(D-86) 


where 


a ri( x i 

‘max 'max y 

a, = 1 / 2 + iXiM/x- -X: i\/4 

‘ \ ‘max ‘max ‘/ 
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D.7 EVALUATION OF THE WAKE INTEGRAL FOR THE OUTER SOLUTION 


We must also include the contributions from the doublet wake integral in the outer solution. 
From equation (110), page 68, of NASA CR-2257, we have 


iA *t f° -'“( x '-v) 

* w = — y e 


x i,+i 


4 >y x dXj 


(D-87) 


iA <p. 
4?K 


cx 

7 


-ICO 


( x_x ii+i) 


x ij + l 


P ydx' 


From the bottom of page 68 this takes the form 


^w 


iXj "VK y 


•/ 


■ico^x'-x^ + j^ -iMXjx' 


iX]Mx 


Hj (Xjr)dx'/rXe (D-88) 


2 2 

where r = y(x - x') + Ky , Since the infinite integral is very slowly convergent, we divide 
it into two parts and change the contour of integration of the infinite integral. Thus 


iMXi x 


iX jYk ye 
Av = - a A ^t 


oo X 

ll 

x x ij + l 


-ico^x'-x^ + j^ -iMXjx' 


H,(x ir )dx' / r >(D-89) 


iXjtKy \ iw ( x i 

Ay> t <e ' | e 

x 




-ico(x'-x) / (3 


Hi(Xi r )dx' /r 


icoxj ,i + iMXi x x . 2 

1+ f -iajx u \ . , , 

+ e / e Hj( Xjrldx / r 

x ii + l 


Now by translating the range of integration of the first integral we obtain 


- X i Ky / e 


CX 

/ 


-ico(x'-x) / (5 


OO 

H l(x 1 r)dx'/r = ^- f e ^ H 0 (xjr)d^ (D ’ 90) 


q/l 2 

where r = y£ + Ky . In reference 1, this integral was made more convergent by changing 
the contour of integration so that we obtain real negative exponentials. The integral then 
takes the form 
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a 


(D-91) 


ay Jq 


f e" 10 ^ H Q (X 1 r)d^ = XjtKy] i / 

J 0 / 0 


IT / 2 2 

-cu’VKy cos 0 / P 


J j^Xjl/Ky sin 0^d0 


7T / 2 


f -wfy c° s 0 / /3 / -— \ 

y o e Y j y vKyXj sin 0j 


2 \ -oj'iK y cosh 0 / (3^ 


7T 


K| ^Xjl/Ky sinh dj 


dd 


2 /* -cdVKy cosh 0 / )3^ 


V 


7T / 2 


Kj^Xjl^Ky sinh 0^d0 


which is computed numerically for each value of y. 


The remaining integral is over a finite range and becomes 


l 


ioox' / j3 


x ii + l 


H 1 (X 1 r)dx/r=-y f 


icox' / ^ 


£j(u)dx' 


where 


x ij+l 


x i r ,2 2~| 

u = — [( x - x') + Ky J 

4 


This integral must also be evaluated numerically. 


Finally, combining equations (D-91) and (D-92) into equation (D-89) yields 


(D-92) 

(D-93) 
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Av“ 


iX jyAip 


1 H 


IT / 2 


-co^Ky cos 6 / i 6 


J^XjYKy sin 0^d0 




p / 2 


-cdVKy cos Oil 3 


Yj^Xjl/Ky sin dj 


(2\ -coVKy cos 0 / j3 / \ 

+ Me K^Xf^Ky sinh 0 j 


d0 


(^j J e'° J ^ y C ° Sh ^ ^ Kj (x jYKy sinh 0 ) d0 j 
' it / 2 


x • , , 2 
-lojx / 13 


iX] 2 YKy iwx ii+1 +iMXjx /* e‘ "' r KjOOdx' 

T M e 


x ij+l 


or 


where 


= A ^t G ( x >y) 


^ 1 Y ( ! -coVKy cos 0 / j3 , <- 

G(x,y) = - — — Si J e J J (X 1 VKy sm 0 d0 




7T / 2 


-toYKy cos 0 / jS 


YjCxVKy sin 0) 


2 -coYfCy cosh 0 / j3 / n/ — \ 

— e K^Xjl/Ky sinh 0 j 


d0 


2 /*°° -coVKy cos 0 / 13 


A 

iX^VKy i^x ii + 1 +iMXjX /* 
e J y • 


K j ^X jYKy sinh 0^ d0 ] 

2 


-itox' / 13 


x i, + l 


£j(u)dx 


and 


A 1 f 2 2 ’ 

u= -~*|(x-x') + Ky 


(D-94) 


(D-95) 


(D-96) 
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Since we have assumed steady symmetric flow, the unsteady perturbation potential satisfies 


^ij m +l " 


m 


^m +2 ^m" 1 


then on the airfoil 

<Vi= -^ii m + Si(^j m -l-0u m ) 

- c s2 (*>a m - - ( d sI F i U + d s2 F i L ) 

= ■ ( c sl + c s2 + 2 )*>ij m + ( c sl + c s2)vij m -l - ( d sl F i U + d s2 F i L ) 

At the trailing edge the Kutta condition requires 

Ai/?j x + iccAi^j = 0 

at x = Xj j , from equation (37) of reference 1 , we have 

c li(A^i 1+ l - &Pi x ) + dii 1 (A^ i - A^.j) + icoA^j = 0 
Solving for A^ + j yields 

Ay> t - A^ ii+1 = A<P- h (l - d Hl / Clil - iw / cj^ + ^d^ / Clij)*^-! 
Using equation (D-98) to define A xp- and A^j j yields 

M = [l - (d 1 . i + iw)/ Ci,] [- (c„ + c s2 + Jm 
+ ( c sl + c s2)*>i,,j n| -l - ( d sl F i, U + d s2 F i, L ) 

+ ( d,i l /Cli l)[" (° sl + c s2 + 2 )*i r l j m 
+ ( c sl +c s2)^i,-l,j m -l - ( d sl F i r l U + d s2 F i|-l L )_ 

Hence we write 


A ^t - + h 2%j m -l + h 3^ r lJ m + Vij-IJm-l + R 


(D-97) 

(D-98) 

(D-99) 

(D-100) 


(D-101) 
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where 

h l = - ( c sl +c s2 + 2 )( 1 - d lij / c li! - icj / c lij) 

, / \ / , . , x (D-102) 

h 2 ~ ( c sl + c s2) ( l - d li / c lij - 1C0 / c lij) 

h 3 = ‘ ( c sl +c s2 + 2 )( d li 1 / c lij) 

h 4 = (^ c sl + c s2) ( d lij / c lij) 

R = - ( ] ~ d lij / c 1 i j ' iw / c lij) (dslFij 17 + d s2 F ij L ) 

- ( d li / c li)( d sl F i -1 U + d s2 F i -1 L ) 

Substituting equation (D-101) into equation (D-94) yields the following expressions for the 
induced flow from the wake at the point xj, yj : 

^w = H lij^i 1 J m + H 2ij^i 1 J m -l +H 3ij^i 1 ,j m -1 +H 4ij^i r U m -l (D-103) 

where 

H 1 ij = h l G ( x i ,y j) (D-104) 


and similarly for the other quantities. 

D.8 DERIVATION OF THE COEFFICIENTS FOR 
A MESH WITH TWO AXES OF SYMMETRY 

To reduce the number of integrals to be evaluated, we consider two lines of symmetry for 
the rectangular mesh region. We define 

^d( x n - x i’Yj) = ^dn^j) = ^dnij 

M X i’ y n- y j) = M X i’ y j) = ^nii 

with similar relations for \j/ un and ^ rn . For i max even, we see from figure 7 that the follow- 
ing relations hold for points on the left and right boundaries: 

^dnij ^di max -n+l, i max ,j 
^£nlj ^ rn i m axj 
^rnlj = ^«ni max j 


Similarly, 
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Figure 7— Illustration of Equal Values of the Integral for Points Associated 
With the Upstream and Downstream Boundaries for a Grid 
With Two Axes of Symmetry 


* 


^dn2j ^di max -n+l, i max -l ,j 
^£n2j ~ ^rni max -l,j 


^rn2j ^£ni max -l,j 

On the lower boundary we have from figure 8 

^ dni 1 ^ di max -n+ 1 , i max -i+ 1 , 1 

^£nil _ ^rni max -i+l,l 

^mil _ ^£ni max -i+l,l 

^dni2 ^di max -n+l, i max -i+l,2 

^ J2ni2 ~ ^ rni max -i +1,2 

^mi2 - ^fini max -i+l ,2 

Because of these equalities we need to compute only the left boundary integrals and the 
lower boundary integrals. For the left boundary, we have 

^dnij ~ ^di max -n+l,i max -i+l,j 1 ~ ^ an( ^ j " 2,3,..., j m and n - 2,3,..., i max -l 

^ finij _ ^rni max -i+l,j for i - 1,2 and j - 2,3,..., j m and n = 2,3,..., j m 

^rnij “ ^fini max -i+l,j for i = 1,2 and j,n = 2,3,..., jm 
On the lower boundary we have 

^dnij = ^di max -n+l, i max ~i+l, j for j = 1,2 and i = 2,3,..., i max / 2 and n = 2,3,..., i max - 
^£nij “ ^rni max -i+l,j 

^mij = ^eni max -i+l,j for j = 1,2 and i = 2,3,..., i max / 2 and for n = 2,3,..., j m 
The total number of \p integrals to be evaluated are: 


1. On left boundary 



O N 

o 



Figure 8 —Illustration of Equal Values of the Integral for Points 
Associated With the Lower Boundary for a Grid 
With Two Axes of Symmetry 


X 


V 


^dn : 2 (jm ^(‘max 2 ) 

^£n : “(jm " (-*m " 

^rn : -(j m “ ') (im ~ 

Total = 2(j m - l) [~i max + 2(j m - l) - 2] 

2. On lower boundary 

^dn • 2 ('max I 2 ~ 0(‘max ~ 2 )~ Omax ” 2 ) 

^Cn ' “(*m ” ')(‘max I 2 “ ') ~ (^m " (‘max - 2 ) 

^rn • (irn * (*max " 2 ) 

Total on lower boundary = (i max - 2^ j^2 ^j m - l) + i max - 2 J 
Combining the two totals yields the following total number Nj of integrals to be evaluated: 

N I = [‘max + 2 (jm' 0 ’ 2 ~f 

D.9 DERIVATION OF THE MATRIX ELEMENTS OF THE 
SYSTEM OF EQUATIONS TO BE SOLVED 

For the sake of completeness we write down the difference equations whose coefficients of 
the ipjj form the elements of part of the matrix. The present program with simpler far-field 
boundary conditions can be coded to compute these coefficients with small modification. 
At elliptic points we have, in the notation of reference 1 , 

- (aj + bj + Ej + E 2 - Qjj^jj + bji£jj +1 + E^ i+1 j + = 0 

for j = 2,3,..., j m and i = 2,3,..., i max -l 

At hyperbolic points, 

a j^ij-l - ( a j + b j - % - Q ij)^ij + Vij+ 1 ' ( E 3 + E 4) *i-l J + E 4^i-2,j = 0 
for j = 2,3,..., j m and i = 2,3,..., i max -l 
For j = j m , we have = - <^y , and the two equations become 

a j^ij-i - (aj + 2bj + Ej + E 2 - Qy) ^y + E ]( p i+1 j + j = 0 

¥>ij-l ‘ ( a j + 2b j ' E 3 - Q ij)^ij - ( E 3 + E 4)^i-1 J + E 4^i-2,j = 0 
For j = j m , we also have boundary conditions for ig < i < ij and jump conditions for i > i j . 
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The coding is already written for the difference equations and we need only consider the 
boundary conditions on the mesh edges. Thus we add to the difference equations the 
boundary conditions 

J m / + iXjMxj + iXjMx 2 \ 

iPi; = ' 5 '' <*n e + \^£ n 2j e ) 


^lj " E ro y J2n 1 j e + 'Pn 

n=2 L 

iXjMxj iXjMaj 
~^£nlj e ” ^nj e 


nj e (x-x 2 )/(4VK) o gn 


'max - ^ 


+ E ^dnlj e ‘ a dn + Y ^rnlj e ‘ * % + ^wlj 

n=2 n=2 

j ^,3 ,..., 


iXjMxj iXjMxj 


+ iXjMxj + iXjMx2' 


^finlj e + ^2n2j e 


iXjMx 2 iX^Ma^ 

'^n2j e ~ + 5 nj e (xj -x 2 ) /(4 Yk) a £n 

'max -1 iXjMx 2 

+ E ^ dn 2j e ff dn 
n=2 

J m iXjMx 2 

+ E ^rn2j e a rn + ^w2j 

n=2 


for j = 2,3,..., j m 

Here the function i p without the plus or minus subscripts designate the combination \ p 

$ ' t ■ 

iXjMxj ( -^m Jm 

^il =e )E ’/ / g n il a gn + E ^rnil^rn 

n=2 n=2 

^ax - ^ n , , 


■E Un« + C a )n-*ann 

n=2 

^ 8 ni(y 1 - V2 ) / 4 ] a dn | + ^Wil 


for i 2,3,..., i max -l 
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^i2 = e 


iX|Mx| ( J m 


J 


m 


Z ^£ni2°Kn + Z ^rni2 a rn 
n=2 n=2 


'max - ^ 


E K 

n=2 


^dnil + ^dni2 ) / 2 " ^dni2 


^ 5 ni(yi "V2) /4 ]°dn| + *wi2 


for i = 3,4,-, i max -2 


iXjMxj 


max 


J 


in 


*max ^ 


j e 
‘max’ 


£ <//£ni maxJ afin+ ^ ^ dni maxJ adn 
n=2 n=2 


Jr 


in 

+ z 

n=2 


“1 ’Z' 


+ 

rni maxj 


iXi Mx 


'max 


iXiMx 


+ ’/' rni -1 j ® 
max ,J 


'max - ' 


iXiMx: iXiMa^ 

1 ‘max „ 1 L 

’/'rni j e ^ ^nj 

‘max’ J 


( x i " x i -]) 
\ ‘max ‘max / 


/(4Vk) 


°rn ^wi j 
111 wl max J 


for j = 2,3,..., j 


m 


$1 - 1 i ® 

‘max 1 > J 


iXiMx, _i ( j 
‘ ‘max ) 


m 


*max - ^ 


\h 

s n=q 


Z ^Kn, i max -l, j°Cn + Z ’/ / dni Triay -l, j a dn 


n=2 


Jm 

+ z 

n=2 


iXjMxj 


max 


iXjMxj 


fl V^-maxJ 


max + ' 'max 

J~ ’/'mi “1 i e 

1Iu max ‘> J 


’Z' rjl i -1 j ^ 

““max 1>J 


iXtMX; t 
1 ‘max 


iXjMa2 

- 5„ ; e ^x ; - x 


nj 

j 2,3,..., j jyj 


J max x max 


/(4VK) + ^ wi .j j 

J W1 max ,J 


Here 


also. 


«q = 1 / 2 + iMXj (xj - X 2 ) / 4 

a i = 1 /2 + iMXi^Xj -x- A /4 
‘ 1 y 'max ‘max J 


We notice that there are two expressions for and -12' ^ ie ^ wo ec l ua ti° ns must be 


max 


equal; hence we obtain for ^2 
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j m r + 
£ L^O^Cn 


£ 

n=2 


iXjMxj 


12 e + ( G '0 " 0 ^£ n 22 e 


iXjMx2 


+ S n2 e 


iXjMa j 


( x l -x 2 ) /(4YK)] 


a Cn 


+ e 


iXjMx 2 l max”^ / + 


nax - / + + x 

£ V^dn22‘ ^dn2V / 2 


11=2 


-^K 6 n2(yi-y2) 


and for .m 

^ m o v * ^ 


ff dn = 0 


max 

J m 


£ 

n=2 


iXiMxj 


_ + 1 ‘max /_ \ + 

G 1 ^rn i 9 e + ( G i " 1 j'/'mi -1 2 e 

In > ‘max," \ 1 / IIU max 


iXlMX: _1 

1 hnax 


" 5 n2 e 


‘max, 
iXjMa2 


fa -xj A/(4Mk) 

\ ‘max ‘max ) 


'm 


iXiMx 


^tn nv 1 


+ e 


^W 1 ‘ max 


n=2 

' S "Wx-l W ( y ‘- y 2) /4 


£ I rdni max -l, 2 ^dni max -l, 1 ) / 2 


°dn = 0 


D.10 THE DEFINITION OF THE BASIC VARIABLES AND 
THE FORMULATION OF THE MATRIX 


We write for the equations ^ 

yi a nm\ ~ ^nr — ^> 2 >-"> N 
n=l 

Let x n - ^ or n “ (i " 1 )Jm " ^ + J f° r x < ' hnax an< ^ n ~ ( 4nax “ 0 fa ” 2 + ^ or * “ faax - 
Let Np be the total number of potential variables. Then 

Np “ faaxfa “ 2 

The total number of ag n variables is j m - 1 (see fig. 8) 

The total number of a rn variables is j m - 1 
The total number of a ( j n variables is i max - 2 

Hence the total number of variables o n is 

2 (fa ” 0 + faax ” 2 

Combining this total with the total number of potential variables, we obtain for the total 
number of variables N 


W 


% 
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X ^max^rn + 2 0m ” 0 + *max ” ^ 


To summarize the previous discussion, we identify the variables x n as 


X n = ^ii for n = (i - l)j m + j - 1 for i < i 


max 


and 


n (*max " l)-im + j " 2 * or * ^max 
X n = ffgk forn = N p + k-l, k = 2 , 3 ,...,jm 

X n = a rk forn = (N p +j m -l) +k-l, k = 2 , 3 ,..., j 


m 


X n = a dk for n = ( N p + 2 J'm - 2 ) + k ” 1 • k = 2 > 3 >- J m 

We now consider the numbering of the equations which make up the matrix system. The 
number of difference equations for the inner solution is 

X d ” Omax ~ 3 )0m " *) 

The number of < p boundary conditions is 

X bc “ ^( J m " * ) + 2 (*max ~ 2 ) ' 2 
The total number of equations is then 

X e - ^max^m + 2 (^m ' 0 + * max ” ^ 

We therefore require two more equations to complete the system. These equations are 

provided by making the relations for 1^99 and p- 1 9 equal. 

^ hnax 1 ^ 

We now define the equation numbers and the corresponding matrix elements. 


1 . Equations numbered m = 1 , 2 ,..., are the difference equations of the inner solution 
and 

N d = (jm ” l)(*max ~ 2 ) 

2 . Equations numbered m = + 1 ,..., + j m - 1 are the boundary conditions on 

J “ 2 , 3 ,..., jjyj 

Nbi = Nd + j m “ 1 is the number of difference equations + the number of boundary 
conditions. 

3 . Equations N^j + 1 to + j m - 1 are the boundary conditions on i/^j , j = 2 , 3 ,..., j m 

N b2 = N bl + Jm - 1 is the nu mber of difference equations + number of <pjj + number of 
</>2j boundary conditions 

4 . Equations m = Nr.9 + 1 to N h 9 + j - 1 are the boundary conditions on 1 ; 

uz oz m imax -1 J 

boundary conditions. 



N^s - + j m -1 is the number of difference equations + the number of p^ + the 

number of ^ 7 : + the number of p-. 1 , boundary conditions. 

J ‘max J 

5. Equations m = N ^3 + 1 to N^ 3 +j m - 1 are boundary conditions on i/ 3 j j, j = 2,3,...j m 

rnsx J 

N ^ 4 = N ^3 + j m - 1 is the number of difference equations + the number of sP j j + the 

number of p->\ + the number of pi 1 + the number of p- ; boundary conditions. 

‘max ‘max J 

6 . Equations m = N b4 + 1 to N b4 + i max - 2 are boundary conditions on 

i - 2,3,..., i max - i 

^b 5 ~ ^b 4 + Jm ~ i i s ^ ie num ber of difference equations + the number of + the 

number of ^ 7 : + the number of p- 1 : + the number of p- + the number of 

“ J 'max J ‘max J 

boundary conditions 

7. Equations m = + 1 to + i max - 4 are boundary conditions on 

^i2> * ~ 3,4,..., i max - 2 

The total number of equations is 

^ ~ Nb5 + *max ^b4 + ^max ~ ^ 

N = N d + 4(j m - 1) + 2i max - 6 

_ ‘max J m + 2 0 m ~ 0 + *max ' ^ 

We require two more equations, since the total number of variables is 

‘max^m + ^(hn " i) + ‘max “ 4- These equations are obtained from equating the two 

relations which give p->n and also give p\ 17 . 

‘max ' 1 z 

Since the wake integral is involved in all equations greater than m = N 4 , we need to identify 
the p variables associated with it; we have 

^wij = * ( H 1 ij^i 1 J m + ^ijjm-l + H 3ij^i r l, j m + H 4ij^i r l, j m -l + R ij) 

Let nj = (ij - l)j m + j m - 1 = ijjm “ 1 

and n 2 = (ij - 2 )j m +j m - 1 =(ij - l)j m - 1 
then 

^wij - “(^lij-^nj + ^2ij^nj-l + ^3ij^n 2 + ^3ij^n 2 -l^ 
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D.l 1 FORMULAS FOR MATRIX COEFFICIENTS 
REQUIRED BY OUTER SOLUTION 

We shall now write the equations for the various coefficients of the matrix a nm . For the 
equations m = 1 to N d , the coefficients are for the inner finite difference solutions and are 
described in reference (1). We now formulate the boundary condition coefficients a., m 
resulting from the source distribution of the outer solution. 


1 . Equations m = N d + 1 to N d + j m - 1 ; 

ipjj boundary condition 

m = N d +j-l, for j = 2,3,..., j m 


a. Coefficients of ipjj 

a M,N d +j-l = 1 


a ni ,m = - H llj 

a n r l,m H 21j 

a n2>m _ ~ ^31j 

a n2*l,m - " ^41j 


when n j = (i ] - ~)(j m - 1 ) + jm and n2 = ( ij - 1) (jm - 1) are the variables associated with 
the potentials about the trailing edge. 

b. Coefficients of variables, k = 2,3,..., j m 

/ + iXjMxj + iXjMx2 

a N p +k-l,m = a 0\^£klj e + ^£k2j e 

iXjMxj iXjMaj 
- ^Cklj e - 5 kj e ( Xi -x 2 )/(4Vk) 

c. Coefficients of variables, k = 2,3,..., j m 

iX^Mxj 

a N+k-l,m = Vi'rklj e 

d. Coefficients of variables, k = 2,3,..., j m iXjMxj 

a N r +k-l,m - ^dklj e 

e. Right-hand side R m = R j j 

+ 

Here the functions are defined in equations (D-16) and (D-17) and \ jj without the 
± signs in this section is understood to be the combination \p + - \jj~ . 

2. Equations m - N^ j + i to N b j + j m - 1 ;^ 2 j boundary conditions 

m = N bl +j-l, j = 2,3,..., j m 
a. Coefficients of ipy variables 
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a j m +j-l,m 1 


a nj,j ' H 12j 


= - Hoo; 


nj-l,m ” 22j 


a n 2 ,j “ H 32j a n 2 -l,m"- H 42j 

b. Coefficients of a fik variables / + iX]Mxj + iXjMx2 


for k = 2,3,..., j m 

c. Coefficient of a rk variables 

for k + 2,3,..., j m 

d. Coefficients of variables 


a Np+k-l ,m “ a 0 \^fik 1 j e * + ^Ck2j e 

iXjMx 2 iXjMaj 

* ^k2j e + 5 kj e ( x 1 - x 2) /( 41 ^) 


a Nn+k-l,m ” ^k2j 


a N r +k-l,m _ ^dk2j 


for k = 2,3,..., j m 
e. Right-hand side R m = R 2 j 

3. Equations m = N^ 2 + 1 to N^ 2 +j m - 1 ; _j j boundary condition 


'max 

m = N b 2 + .i- 1 for j = 2,3,..., j m 


a. Coefficients of y?jj variables 


3 imax' 2 Jm+H,™ 


= 1 


= - H 




11 r 1 ,m 2 ^max”^’j 


n l’ m “max 

a n 2 ,m " ' H 3i max -U a n 2 -l,m = ~ H 4i max -l,j 


A 3i -1 i 
m aY 5j 

b. Coefficients of crg k variables 


a N n +k-l,m-^Cki max -l,j e 


iXiMxj _i 

4 ov 4 


max 


for k = 2,3,..., j m 
c. Coefficients of o rk variables 
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a N c +k-l,m- a l\^rki max j e 


iX. Mx: 

+ 1 ‘max 


iXiMx; i 
max + 1 ‘max ‘ 

+ lKb-j _i e 
1K1 max 1 


* rki max- 1 J 6 


+ 5i,; e ( 


max 'max 


ax"')' 


for k = 2,3,..., j 


'm 


d. Coefficients of variables 

a N r +k-l,m _ ^dki max -l,j 

e. Right-hand side R tn = R ; 1 : 

‘max 


iXiMx= i 
1 ‘max 


for k = 2,3,..., i max - 1 


4. Equations m = N b3 + 1 to N b3 + j m - 1 ; ^ j boundary conditions 

max 

m = N b3 + j - 1 for j = 2,3,..., j m 

a. Coefficients of <pjj variables 

'max"' j m + j'^ ,rn 

n l ,m "max’j nj-l,m 2i max ,j 


* n 2 ’ m 3 i max>j V 1 -™ 4 'max’J 


b. Coefficients of variables 


iXjMx. 


a Np+k-l ,m “ e 


^CkimaxJ 


for k = 2,3,..., j. 


c. Coefficients of variables 


iXjMxj 


I + max 

a Ng+k-l ,m - a l V^rki flY j e + 7ki _„-l , j e 


iXiMx: i 

max 


iXjMx- iX Ma 2 

^rki i e maX + §ki e -x: A/(4VK) 

rkl max J Kj ^ 'max 'max ‘J v 


for k = 2,3,..., j, 


d. Coefficients of variables „ 

1 Ifns 


1 max 

a N r +k-l,m “ e ^ dki r 


for k = 2,3,..., i r 
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e. Right-hand side R m = R: .- 

111 'max J 


5. Equations m = N b4 + 1 to N b4 + i max - 2 - 

m = N b4 + i - 1 , i = 2,3,-, i max -l 

a. Coefficients of variable 

a i- 1 i .m — ^ 

* i Jnr 


a n j ,m ~ ^ 1 i 1 


a n-),m *^3il 


b. Coefficients of agj, variables 
for k = 2,3,..., j m 


a Np+k-l,m e 


a ni-l,in ~^2il 


a n2-l,m ~^4il 
^Ckil 


iXjMxj 


c. Coefficients of variables 

a N c +k-l,m = e 


iXjMxj 


^rkil 


k = 2,3,..., j 


m 


d. Coefficients of oy k variables 

a N r +k-l,m = e 1 [(^dkil + ^dki2)/ 2 -'/ / dkir'^ 5 ki( y l _ y2) /4 ] 
tor k — 2,3,..., ijyj a ^ ” 1 

e. Right-hand side R m = Rjj 

6. Equations m = N b g + 1 to + i max - 4 

m = Nb5 + i-2, i = 3,..., i max - 2 
a. Coefficients of v?y variable 


a i-1 j m +l,m 


= 1 


l nj,m "^li2 

a nj-l,m ' ^2i2 

! n2,m ~ ” ^3i2 

a n2”l ,m _ “ ^4i2 


b. Coefficients of variables 


70 



iXlMX; 


a N p +k-l,m ~ e ^Cki2 

for k = 2,3,..., j m 

c. Coefficients of variables 

iXjMx- 

a Ng+k-l,m = e ^rki2 

k 2,3,..., j 

d. Coefficients of variables 

iXjMxj r/ + + \ 

a N r +k-l,m = e [_V dkil + ^dki2 )l 2 “ ^dki2 + *^ 5 ki 1 ~^2 ) / 4 

for k 2,3,..., ^ 

e. Right-hand side R m = Rj2 

7. Equation m = + i max - 3 

Matching of two relations giving ^2 

a. Coefficients of ^ variables are zero. 

b. Coefficients of crg^, k = 2,3,..., j m 

+ iXjMxj + iXjMx2 
a N p +k-l ,m = % ^£kl2 e + ^fik2,2 e 

iXiMaj 

+ 5 k2 e (x 1 -x 2 )/(4Yk) 

c. Coefficients of variables are all zero. 

d. Coefficients of a^, k = 2,3,..., i max - 1 

iXjMx 2 ry + + \ -i 

a N r +k-l,m = e Lvdk22-^dk21 / )/ 2 ' ”^ 5 k2 ( y l ' y 2) / 4 J 

e. Right-hand side R m = 0 

8. Equation m - N b5 + i max - 2 

Matching of two relations giving i 9 

^max 

a. Coefficients of i/>y variables are all zero. 

b. Coefficients of crg^ variables are all zero. 
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c. Coefficients of a^, k = 2,3,..., j m 

a Nn+k-l,m ~ 0i l^rki Triax 2 e 


iXiMx: 

1 l , 


iXi Mx. 


+ (“l-')*rki max -1.2'= 


l MX i -1 
1 'max 


+ S k2 e 


iXjMa2 


i " x i -1 
‘max 'max 


d. Coefficients of o^, k = 2,3,..., i max - 1 


) /(4Vk) 


iXiMx.- i 
1 'max ' 


1 'max -1 )/ + \ 

a N r +k-l,m * « |^dki max -l,2 - *dki max -l,l) 

' 8k imax-'^( y l- y 2) /4 } 


e. Right-hand side = 0. 


The integrals i// ( j n , \//g n , and \p m can be expressed in the form of a single function of 
several variables resulting in considerable saving in coding. These integrals are represented by 


* d ” (X ’ y) = “4 


-j f n ( x ' -x n-l) C 0( u d) dx '/ 5 n 

n [ J* n -\ 


/* x n+ 1 

/ ( x n+l -x ') C o( u d) dx '/ 5 n+l 

-/x„ 


+ i 2 

^gn( x >y)" - “ ~ ~ 

4i 5 n+1 + 5 ] 


- U V-Vn- 

n U y n -l 


l) { o( u i) d y'^n 


+ J v (y n +i -y') fi o( u e) d y'/ 5 n+i 

jy n 


^rn( x -y) = - 7: ~ 

4i 5 n+1 +5 


-{ f n (y'-y n -i) fi o( u r) d y'/^ n 

n U y n -l 


/* y n+ 1 

J v (yn+i - y') e o( u r) d y'/5 n+ i 

j y n 
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where 


x n _x n-l’ yn'^n-1’ 


u d = V 


2 2 ] 

(x' - x) + K(y + b) J / 4 


u g = M 2 [( x ' a i ) 2 + K (y ' y') 2 _ / 4 
2 T 2 2 

= Xj K[(y-y') + (x-aj) / Kj / 4 

u r =x i 2 [(x- a 2 ) 2 + K (y -y ') 2 / 4 
= Xj 2 K [(y - y') 2 + (x - a 2 ) 2 / K 


( 2 ), 


(D-108) 


(D-109) 


(D-110) 


(D-lll) 


and Cg(u) = Hq (f) with f = 2 Vu . The functions i//g n , and \fr m can be written 
down by replacing y by -y. 


Let 


u = ju 


2 2 ’ 

(x'-p) +£ 


(D-112) 


and define the function 


{/ "(x'-x n .,)%(lOclx'/S n 
x n-l 


+ 



- x')£g(u)dx' / 5 n +l 


(D-113) 


where 5 n = x n -x n _j as before. 

By comparing u in equation (D-l 12) with u^, Ug, and uj. in equations (D-109), (D-l 10), and 

(D-l 1 1) and comparing equations (D-105), (D-106), and (D-109) with equation (D-l 13), we 
see that 


tf'dn = ^( X l 2 / 4 > VK(b + y),x,x n ) 

(D-l 14) 

^£n = ^( X l K/4, ( x -aj) / VK,y,y n ) 

(D-l 15) 

*rn =4*l 2 K/4, (x-a 2 )/YK,y,y n ) 

(D-l 16) 


Hence the subscripted quantities become 
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4nij = v£/ (Xj 2 / 4, VK(b + y j )) x i ,x n ) (D-l 17) 

^Cnij = ^(x^K/4, (xj -a) / VK,yj,y n ) (D-l 18) 

^rnij =^(^l 2 K/4, (x r a 2 ) /l^, yj ,y n ) (D-l 19) 

^dnij =^(xi 2 /4,YK(b-y j ),x i ,x n ) (D-l 20) 

^nij “ ^ (x 1 2 K / 4, (xj-aj) /VK,- yj ,y n ) (D-121) 

^mij - * (x^K / 4, ( Xj - a 2 ) / VK, - yj ,y n ) (D-l 22) 
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APPENDIX E 


AN OBLIQUE COORDINATE SYSTEM FOR 
SWEPT AND TAPERED WINGS 


Consider a vector function F with x,y,z, components F j E2E3 ; then the divergence of F 
under the transformation of equation (G.l) becomes 

- 9 F 

V-F= — •V6+F 2n + F 3r| (E-l) 

where g = (g x ,g y ), and F j JF2T3 are the x,y,z components of F. Expressing the opera- 
tor in conservation form yields 

9 ' ' (E- 2 ) 


V-F= (F -Vg) + F 2 „ + F 3b - F • 77 (Vg) 


Now 




Vg = (l / c(y), £c' / c - x' £e (y) / c ) 


Here c may be written as c(r?) and we find that 


We also have 


9 £ 


Vg = ( 0 , - c' / c) 


F= ^uv? lx - 2 ico»p 1 I e ,y> ly 


or 


substituting the transformation yields 


(E- 3 ) 


(E- 4 ) 


(E- 5 ) 


(E-6) 


Since the linear transonic small perturbation equation for unsteady flow can be written in 
theform VF+q ^=0 


we obtain from equations (E- 3 ) through (E-6) 

3 0 

^[Sx(^l { 8 x -2i^, / e ) + *y(*l„ + *y*l t )] + ^(% + Vl f ) < E 


;-?) 


The first derivative terms must 
Thus 


+ (c' / c)^?^ + gy^]^) + q^i 


= 0 


\ u <> / 

be changed to reduce the equation to conservation form, 
c' 9 /c' 

T 


r\ 817 \c 


*1“ 7 


g y% 3^ (V^) 
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The conservation form of the differential equation then becomes 


a_ 
at 

G 

+ — 


(g x “ u + gy 2 )<Pi + (gy + 2ico / e)«pj + g y ^! 


(E-8) 


3 ” 

-|0+c'/c)*> 1)i + g y *> 1{ J+¥>, 


+ [q + c' / c - (c' / c)'] j = 0 


rr 


For the coefficient of v?j we see that the second derivative of the chord in the span wise 
variable must be continuous for the conservation form to be valid. If the nonconservation 
form is used, there is not this restriction on the planform and we obtain 


g x (a / at) jsx^ig - 2kj ^i / e J + g y ( 9 / a ^)(^i T? + 8 y^l^) 

+ (a/3rj)^)j^ + gy(Pi^ + (Pi^ + q«Pi =0 

The condition that the equation be hyperbolic for both forms is 

2 2 

g x U + gy <0 


(E-9) 


(E-10) 


The root chord of the wing must be a plane of symmetry and we must impose the conditions 

(E-l 1) 


that ipi = 0 at y = p = 0. In terms of the £,r? variables this becomes 

y 


<P 1 + gytf>ig = 0 for y = 77-0 
Since this term is zero we must also have, for small p 2 , 

(^r? 8yiPl f)r? = -t?2 /2 (^t? gyV?1 t)r? = r? 2 /2 

and the difference form of the rj derivative becomes, at 17 = 0, 

_ (% + g y^)?? = t? 2 /2 


d 


(% + g y%) 


(E-12) 


' *j/ r? 2 

If we introduce the quantities 

M = g x = 1 / c(r?) , ^ = g y = - £c'(r?) / c(i?) - xg e (r?) / c(r?) 

then the differenential equation (eq. E-9) becomes 


M(3 / 3$) mu</ 3|^ - 2ico^j / e + v(d / 3$) 


(E-l 3) 


+ 


(3 / 3i?)^i^ + + q^i = 0 


♦ 


* 
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Using the equations (E-l 1), (E-12), and (E13) yields the following differential equation for 
points along the root chord rj = 0 of the wing 

p(8/9£)pu^-2 ku^/e + (pi^ + JVi^ /2 IV2 + 'P\tf + Wl =0 (E-l 4) 

Equations (E-l 3) and (E-l 4) may be differenced in the same way as the differential equa- 
tions in reference 2. Thus for the first derivative in £, we have 

= c li (*i+l, jk “ ^ijk) + d li(^ijk - *i-l, jk) (E-l 5) 

where from equation (H-20) on page 40 of reference 2 

c li = (^i " ^i-l) I [(^i+1 ' ^i-l) (^i+1 ' *i)] (E-16) 

d li ~ ($i+l ' ^i) I [($i+l ' ^i-l) (^i ” *i-l)] 


Similarly 

% = c lj (^ij+1 , k - ^ijk) + d lj(*jjk - ^ij-1 , k) 
where cjj and dy have the same form as equation (E-16) but with rj replacing £. 

Since p is a function of y and v is a function of x and y, we may write 
M( T ?j) = Pj and v at the point (i,j,k) 

From equations (19) and (20) of reference 1 we see that 


2 r 

- Pj [2c i u i+1/ / 2 ji c ^i+ijk -v?ijk) 

■ 2d i u i-l/2,jk(^ijk-^i-ljk) 


2pj (i w / e) Jc j j («p i+ j jk - y> ijk ) + d H (V ijk -^_ 1; jk )J 


(E-16) 


(E-l 7) 


(E-l 8) 


"«[ Cll (% i+ i >jk %jj k ) +dll C\jk %_!,*)_ 


(E-l 9) 
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Applying the formula in equation (E-15) to each term of equation (E-19) we have, finally, 


"5$ % ij k = y ij j c li [ c l.i (^i+1 , j+1, k " ^i+l, jk) + d lj (^i+1, jk ~ ^i+1, j-1, k)] (E . 20) 

+ ( d 1 i ’ c 1 i ) [ c 1 j (^ij+ 1 , k ' ^ijk) + d 1 j (^ijk ' ^ij- 1 , k)] 

- d li [ c lj (^i-1 , j+1 , k - *>i-l , jk) + d lj (W-l , jk - *i-l , j-1 , k)] | 

Similarly, from equations ( 1 9) and (20) of page 40 of reference 1 we obtain 

” T( (* V U) =,, ij [ 2c i"i+l/2, j^i+l, jk "^iik) - 2d i I 'i-l/2,j(* , iik -VH, jk)] (E-21) 
The remaining second derivative terms take the form 

n = 2a j (^ij+ 1 , k - ^ijk) - 2b j (v’ijk ' ^ij- 1 , k) (E-22) 


(E-22) 


^1 ^ ~ 2a k(^ijk+l " ^ijk) ' 2b k (^ijk " ^ijk-l) (E-23) 

The boundary conditions on the wing take the same form as for the cartesian coordinate, 
since the is essentially unchanged from the unswept case derived in reference 1. 

Consider the term in equation (E-14) 

(V‘ Vl (&24) 

We need to express it in terms of the quantities at the grid points. Now 

lpl T ? |r ? = r ?2 /2 = ( lPi2k ' V,ilk ) /r?2 (E-25) 

Remembering that r/j = 0 for j = 1 and rj = pj for j = 2, we see that 

^Ipjp = t? 2 /2 ~ C li(^i+1, 2, k + < ^i+l, 1, k“^i2, k "^il, k) 


(E-24) 


(E-25) 


(E-26) 


+ d li(^i2, k + ^il,k’^i-l, 2, k'^i-1, 1 , k) / 2 

Substituting equations (E-25) and (E-26) into equation (E-24) leads to 
r 1 2 

K + ■Vlf), = , 2 / 2 J 1 ”2 = (^2. k - Vil, k) / 02 

+ ^i,3/2 [ c li (^i+1, 2, k “ ^i+1, 1, k " ^i2, k “ ^il, k)] / 2t ?2 

The assumption of plane wave boundary conditions on upstream, downstream, upper and 
lower boundaries in the cartesian coordinate system yields 
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i£i + icoMyb / (1 + M) = 0 at x = X: 

1 x 1 ‘max 


n 


X 


+ icoM <Pj / (1 - M) = 0 


at x = x i 


-iXj"^<Pj=0 at y = yj 


max 


<p ] +iXj"VK<Pj=0 aty = yj 

In the new variables, these relations become 

+ icoM<Pj / (1 + M) = 0 

- icoMY’j / (1 - M) = 0 

¥>1 + + iXi Vk^i = 0 

V $ 

+ *Vi ~ iXi VK^ = 0 
i? € 

From the form of the second derivative terms and particularly the cross-derivative terms we 
see that the difference equations associated with an interior point involve the 1 1 points in 
figure 4 in place of the usual 7 points in figure 5. For the x = constant line relaxation solu- 
tion, the matrix is still tridiagonal. For a direct solution, the matrix is still sparse and 
somewhat banded. 


Along the wake , the condition that the vortex sheet not support a load is 

+ icoAi^j = 0 

and is thus unchanged from the cartesian form. 
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APPENDIX F 


ROW RELAXATION FOR THREE-DIMENSIONAL FLOW 


For hyperbolic points, additional fictitious time-dependent terms must be added to make the 
row relaxation procedure converge. Under the assumption that the calculations are swept in 
the direction of increasing j and k, the difference equation (A-2) from reference 2 becomes 

2E 3 (4jk - 4-1 , jk) + 2E 4 (4-2, jk " 4-1 , jk) (F .j } 

+ o ( (n) (s)N i on ( (s) (n - 1} \ 

+ 2a yj ^ jj _ 1 ; k - <p ijk ; - 2b yj ^ ijk - *> ij+1 s k ) 

+ n ( (n) (s M on ( (s) (n -^ 

+ 2a zk\^ijk-l “ ^ijk/ - 2b zkWjk " ^ijk+1/ 

(s) 

+ q ijk^ijk “ 

where E 3 = c i _ 1 u i _ 1/ ^ jk - ico c 2i / e, E 4 = d i _ 1 u i _ 3/2jk - iwd 2i / e, 

and the superscripts n and n - 1 denote the results of the current relaxation sweep and the 
previous one, respectively. The superscript s denotes the quantities for which equations 
(F-l) for all i and for fixed j and k are solved. The subscripted variables a, b and c are 
defined in reference 2 on page 68. 


We now introduce a fictitious time derivative related to the iteration by the relation: 

(n-1) _ (n) J (n)\ 

^ijk ^ijk _2vt \^ijk/ t (F-2) 

Introducing underrelaxation by a factor r yields 

(n) _ ( (s) (n-l)\ (n-1) 

(n _l) ^ijk “ r V^ijk " ^ijk y + ^ijk 

Eliminating i^jjk leads to 

(s)_ (n) /r-l\ / (n)\ 

nik - ^ijk - \y) Ax v’iik K (F ‘ 3) 

by means of equations (F-2) and (F-3), the difference equation can be expressed entirely in 
terms of the n tb iteration for the potential. After dropping the superscript n, we obtain 
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2E 3 
+ 2E4 
+ 2a yj 

" 2b yj 

+ 2a zk 


r- 1 


^ijk'^i-l ,jk"( r ) At (^ijk~^i-l,jk) t 


(F-4) 


r- 1 


*i-2, jk - *i-l , jk - [—) At (*i-2, jk - *i-l , jk) t 


^ij-1, k ' ^ijk 


+ At 


>(^k) t 


^ijk 


r- 1 


At ^ijk t -^ij+l,k + At (^ij+l,k) t 


r - 1 . 

^ijk-1 -^ijk + i^-] At (^ijk) t 


-2b 


zk 


^ijk ' ^ r ^ At (^ijk) t ' ^ijk+1 + At (^ijk+l) t + ^ijk ^ijk ' ^ r ) At (^ijk) t 


= 0 


Replacing the difference terms by their appropriate derivatives in preparation for taking the 
limit as the grid size goes to zero yields 

- 2b y jAt(V|j +1) k - Y3jj k ) t = - 2b y j (yj +1 - yj)At y> yt 


Substituting for b y j from equation (A-3) of reference 2 leads to 


- 2b 


2At 

y j At (v>ij+l ; k - ^ijk) t “ - ^7 fyt 


(F-5) 


where Ay = yj + j - yj_j - A similar y> zt term results from the b k term. The terms result- 
ing from the first two terms of equation (F-4) cancel on taking the limit as Ax, Ay, Az go 
to zero. We now proceed to the limit. Neglecting terms of order one and higher in the 
small increments we obtain the following differential equation from the difference equation 
(F-4) 


U *T> 


2ico 






At 


yy l zz Ay 


:riyt 


Ayj 

Azt 


^lztj 


(F-6) 


- 2At 


*lt 


+ qi^j = 0 


b y j / r - a y j(r - 1 ) / r + b zk / r - a zk (r - 1) / r 
h = [ b yj + b zk - ( a yj + a zk) (r " 1 )] A Vj / r 

[u«,, x -2i( M /,V 1 ] x + , 1 yy + Vl z z -2^-^ lyt+2 | % t + <>3%Wl-° < F - 8 ) 


Let 


then 


(F-7) 


The differential equation finally becomes 
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(F-9) 


ui^Ji + y>i + <^i + ( u Y - 2ico / e) <p-\ 

xx 1 yy *zz v ' J x 

At 

" 2 Ayj (^*yt + ^ jk ^’zt + ^'T) + ^1 = 0 


where 


% 


yj+i-yj-i 

z k+l ~ z k-l 


To obtain a convergent operator we add terms (Jjipj + 02^1^. to the terms in parentheses 
and determine j3 1 and (L* to yield a differential equation with the correct time-like behavior 
for the x coordinate. Let r = a^x + a^y + a^z + t then 


*l x _0! l*l + *l x 


V\ ~ a 2^l +< ^1 

l y z l r l y 


<P\ =( *\ 2 V\ +2a lV ?i +¥>i 

2 xx 1 Vr 1 xt xx 


<Pi ~ a 2 ^1 + 2a 2^1 + ^1 

l yy z a tt z *yr l yy 


(F-10) 


*l z = a 3*l T + *l z 


~ ^1 + + <P t 

zz J *rr zt *zz 


^1 * -a l^l +l P\ 

*Xt 1 l TT l XT 


n yt =a m TT + 


n 


y t 


<P\ ~ a 3^1 +{ P\ 

l ZT J l TT ZT 


Substituting equations (F-10) into the differential equation yields 

/ 2 \ 2 
u(ai ip 1 +2ai^>i +^>i +«i if] + 2 « 9 ij?i +<^i 

\ 1 1 tt 1 1 xt X xx/ ^ 1 tt ^ yT l yy 


2 

+ « 3 <P! TT +2a 3 ^ lzr + *>l zz + ( u x“ 2iw/e)(v> 1 x + “^ 1t ) 
- 2 (At / Ayj) ^i rr + ^i yr + 0i («i^i rr + n xT ) + h*\ T 


(F-l 1) 


+ 0jk(«3^1 rr + ^l zz ) + %^l r ] + ^1 "° 

It is necessary to eliminate the cross derivative terms in time r to reduce the differential 
equation to canonical form. This requires 
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uaj - j3jAt / Ayj = 0 


from which 


(F-12) 


«2 - At / Ayj = 0 

0'3-/3 j k(A t/ A y j) = ° 

aj =jS] (At / Ayj)/ u 
a 2 = ( A t / Ayj) 

a 3 = ^jk( At /A yj ) 

Tire elimination of the <^ T terms yields 

«j (u x - 2ico / e) - 2j32At / Ayj - 2/? 3 At / Ayj = 0 

The quantity j3 3 is given in terms of j3j and |3 3 by substituting (At / Ay^ / u for 0 !j . 
Thus 

^2 = ^1 ( u x ” 2i<X) / e ) / 2u " % (F-13) 

In order for the x variable to be time-like the coefficient of must be positive. This yields 
the following relation 


2 2 2 At r .. 

°T u + «2 + a 3 - 2 ^ [ a 2 + ^l a l + ^jk a 3] >0 


Substituting for the oj terms yields 


J 


/ \ 2 
/At \ 

2 , 2~ 

"Wj/ 

01 /u+l+j3j k 


>0 


(F-14) 


(F-15) 


Since u < 0, we have |3j > (- u) ^1 + ) or 

]/l + |3j k 2 , c>l 

02 = @1 ( u x " 2ico / e ) / 2u " % 

% = - (yj+i - yj-i) h + b k - ( a j + a k)( f - *)] / r 

We have now determined the values of j3j and ^ required to establish a convergent operator. 
The differential equation which now must be differenced is 

2i ( w / e )^i] x + % y + n zz ~ 2 [M xt + hn^\ + = 0 ( p - 19) 


(F-16) 

(F-l 7) 
(F-18) 
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Now 


(n-1) _ (n) 
^ijk “ V>ijk 




Eliminating 


(n) 

^ijk 


yields 


(s) 

^ijk 


(n) 

^ijk 


r - 1 


At 


9_ ( (n)\ 

at / 


(s) (n-1) _ 1 

^ijk ~^ijk - r 



At 


3 ^ijk _ / (s) (n-l)\ 

ai r Wjk-^ijk ) 


(F-20) 


(F-21) 


The resulting difference equation then becomes (see equation (F-l)) 


2E 3 (4jk ' 4-1 , jk) + 2E 4 (4-2, jk " 4-1 jk) 

+ 2a j (4-1 , k ' 4 § k ) ' 2b j (4k ‘ 4+ I?k) 


, 0 ( (n) (s)\ / (s) (n-l)\ 

+ 2a k Wjk-1 ~ ^ijk ) ~ 2b k\^ijk - ^ijk+ 1 ) 


(F-22) 



We now consider the case in which the z or k variable is swept in the direction of decreasing 
k. The only terms that change in equation (F-l) are 


/ (n-1) (s)\ / (s) (n) \ 

2a k V^ijk-1 ” ^ijk/ “ 2E k Y^ijk " ^ijk+1/ 


(F-23) 
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c 


Eliminating < p 


(n- 1 ) 


(s) 

by equation (F-2) and ^ by equation (F-3) yields 

r- r 


2 ul 


^ijk-i - At(^ k -l) t -^ij k 


+ 


At 


(^ijk) t 


2 b. 


^ijk 


r- 1 


At (^ijk) t "^ijk +1 


Now 


2a k At ^ ijk - t = 2a k (z R - z k _,) Aty zt 


From equation (A-3) of reference (2) we have 


2a z kAt JjPjj k ~ ‘r’ij k- 1] t 


_ At 

At = ZZ ** 


z k+l “ z k-l 

Taking the limit as Ax, Ay, At yields for equation (F-24) 

/r - 1\ At 

^zz + 2b zkl r ) ^t " 2 a zk A * / r + ^zt 


(F-24) 


(F-25) 


In place of equation (F-9) we obtain 

U ^ 1 XX + % y + ^ 1 ZZ + ( U X ' " 2 ( At / A Vj) (^l yt - ^l zt + ^l t ) +4^1=0 

(F-2 6 ) 

where $3 is now given by 

% = [ a zk + b yj - ( a yz + b zk) < r “ 1 )] A Vj / r ( F ‘ 27 ) 

instead of equation (F-7). Equation (F-26) differs from equation (F-9) in form only in the 
sign of the < p z ^ term. Thus the third line in equation (F-l 2) becomes 


At , 

“3 = - ^jk / 2 


J 


(F-28) 


Equation (F-l 3) remains unchanged but the sign of <^ 3 /3j k in equation (F-l 4) is changed. 
Substitution of equation (F-28) for <* 3 , however, yields equation (F-l 5) unchanged. Thus 
the only change in j3j and j3 2 is the definition of j3 2 by equation (F-27). 


The correction terms for three-dimensional flow do not differ in form from the terms 
derived in reference 2 for two dimensions. The quantity j3. differs only by the factor 

i + ( A yj 1 Az k f 

from the two-dimensional value. The coefficient of the j3 2 term contains the additional 
term contains the additional term from the z derivative (increasing) 

[ b zk - a zk( r " / r 

analogous to the two-dimensional relation 

[ b yr a y j< r -!)] /r 

(This term is given in error in reference 2 as - a y j(r -1) / r .) 
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